CFD and Experimental Characterization of a Bioreactor: Analysis via Power Curve, Flow Patterns andkLa

Abstract

Mixing operations in biological processes is of utmost importance due to its effect on scaling-up and heat and mass transfer. This paper presents the characterization of a bench-top bioreactor with different impeller configurations, agitation and oxygen transfer rates, using CFD simulations and experimental procedures. Here, it is demonstrated that factors such as the type of impeller and the flow regime can drastically vary the operation as in the preparation of cultures. It was observed that the bioreactor equipped with a Rushton generates a kLa kLa of 0.0056 s⁻¹ for an agitation velocity and airflow rate of 250 RPM and 5 L/min, respectively. It is suitable result for the dissolved oxygen (DO) but requires a considerable amount of power consumption. It is here where the importance of the agitator’s diameter can be observed, since, in the case of the two propeller types studied, lower energy consumption can be achieved with a smaller diameter, as well as a much smaller shear cup 2.376 against 0.723 s⁻¹ by decreasing by 4 cm the standard diameter of an agitated tank (10 cm). Finally, the kLa kLa values obtained for the different configurations are compared with the maximum shear rate values of different cell cultures to highlight the impact of this study and its applicability to different industries that use agitation processes for cell growth.

1. Introduction

Several sectors of industry use mixing processes, for example, in the synthesis of products such as oils, cosmetics, pharmaceuticals, food additives, etc. Notably, in the biotechnology area, mixing is a way of providing a correct microbiological environment in which microorganisms can grow with the best living conditions. However, these processes are difficult to control. Even though it is possible to establish the growth conditions in controlled surroundings, the scaling-up to industrial levels of these processes can be troublesome. This is due to many variables, which are linked to specific agitation parameters such as power number, pumping number and impeller type [1]. Many industries have proposed optimization methods for steady mixing processes. However, some small-scale operations that use stirred tanks are emerging, such as two-phase flows and agitation of fluids with high viscosity or non-Newtonian behavior [2]. In addition, other parameters must be considered, such as ideality and regimen in which the agitation process is being carried. Concerning the latter aspect, the assumption of ideal agitation processes is often used for reaction kinetics calculations, which need to ensure a uniform mixing throughout the reactor at any given time.

On the other hand, controlled kinetics are the cases where it is not possible to assume ideality, only being characterized by the mass transfer, where the agitation velocity and impeller type are mostly irrelevant. Additionally, considering solving partial differential equations for momentum, energy and continuity for an omnidirectional flow increases the complexity of the model. Taking these differential equations into account makes problems more challenging to solve to a point where numerical methods are needed to obtain a concrete solution. Some alternatives include finite differences, finite volumes or finite elements [2].

Another aspect to consider in agitation processes is the flow pattern, categorized into tangential, axial and radial. However, for a cylindrical tank without baffles and with any impeller on its vertical axis, the flow patterns at the farthest distances from the impeller will mostly be tangential. In the areas closest to the vertical axis, the impeller type that is used establishes the flow patterns. The axial flow impellers generate a flow parallel to the axis up- or downwards according to its design; the radial flow impellers create a flow perpendicular to the direction of the rotation axis. Many of the impellers available today fulfill the requirements for controlled processes.

Furthermore, in the case of bioreactors, some of the sensors used in the measurements serve as baffles and reduce vortex formation, but they are not as efficient as flat-walled baffles [2]. Stirred methods are widely used in processes where it is necessary to enhance the mass transfer; the bulk generated by the motion facilitates faster and more efficient mixing. Although there is molecular diffusion, also called conductive diffusion, global mass transfer is mainly affected by convective diffusion. Mass transfer is of great importance for bioprocesses due to the oxygen requirement in aerobic cultures. The microorganisms present in the culture consume the available oxygen from the liquid phase, and, in turn, the liquid phase becomes saturated with oxygen from the bubbles, which enter the system through the diffuser. A differential concentration propitiates the mass transfer from the bubble’s interface to the areas where no oxygen is present in the culture [3].

Mass transfer in bioprocesses is affected by several factors that can be classified into two types, biological and operational. The former considers the type of cell required and the substrate used, while the latter considers aspects such as mechanical power, sort of impeller, etc., [4]. This study focused on the operational conditions that ensure the most significant mass transfer from the gas bubbles to the liquid. This was done using Computational Fluid Dynamics (CFD), corroborated by experimental approaches.

2. Literature Review 2.1. Previous Work

Studies to date show the importance of analyzing multi-impeller agitation systems to determine their operational performance. Gogate and collaborators investigated agitation processes with single and multiple impellers [5]. They found that, regardless of the type, the total power will be approximately the sum of the separate powers. Additionally, they affirmed that the flow patterns in a system with multiple impellers are remarkably similar to those with only one if the distance between them does not exceed the diameter of either impeller. In addition, they provided some correlations for calculating _kLawith a different configuration, for example, in two Rushton impellers. Buffo and colleagues analyzed the shear rate in systems with multiple impellers, for three different combinations, obtaining correlations that connect the_kLa , dissolved oxygen (DO), agitation and rheological constants of the fluid to be mixed, determining this parameter of great importance for the cultivation of microorganisms [6].

2.2. Dimensionless Numbers

The impeller supplies movement to the tank contents, integrating the components of the mixture, increasing homogeneity and stability for a prolonged time. Agitation is considered fast when the agitation velocity exceeds 500 RPM, due to the absence of baffles in the bioreactor, the presence of vortexes at high speeds is inevitable, generating a physical restriction in the system. In general, for cases where slow agitation is needed, a speed reducer is used, which is the case in biological processes [7]. The impeller choice depends on several factors, including mixing time and energy consumption. Some of the criteria for agitation are the type of mixture (i.e., suspension, homogenization, dissolution, diffusion and suspensions), the volume of fluid and the physicochemical properties. Finally, the flow regime is determined. This is relevant since, in the laminar or viscous regime, there is little or no introduction of air to the mixture, in contrast to the turbulent regime where there is high oxygen dissolution, directly affecting the impeller power and growth rate of cells in bioreactors. Therefore, the Reynolds number (Re) is determined by the expression in Equation (1), whereρis the fluid density in kg/m³, N is the angular velocity in RPS, D is the impeller diameter inmandμis the dynamic viscosity in Pa·s.

Similarly, the relationship shown in Equation (2) is used to determine power, where_Np (power number) can be found in the literature and depends on the impeller type. According to several authors, the passage from laminar to the turbulent regime is not immediate [8,9,10]. For the case of an agitated tank without deflectors, the laminar regime can be considered for Re lower than ten and turbulent for Re above 10,000 (intermediate Re numbers are defined as in the Transition regime) [11].

Re=ρN^D2μ

P=_Npρ^{N3 D5}

2.3. Oxygen Diffusion

The most common way to determine the transfer of oxygen in a stirred tank is assuming that the mass diffusion occurs in non-uniform, binary mixtures, where the concentration of any component varies along any coordinate of the bioreactor. The molecules move from the area of the highest concentration to the lowest, with the addition of a motion source, e.g., mechanical agitation (mixing is mostly by convective phenomena). Thus, if an area perpendicular to the direction of diffusion (e.g., y-axis) is considered, Fick’s law (Equation (3)) can be used, which states that the mass flow of component A is proportional to the concentration gradient in that direction [3].

_JA,y=_NAa=−_DAB(d_CAjdy)

whered_CAj/dyis the change in the concentration of component “A” along the “y” direction, also known as a concentration gradient. j can be L or G depending on the liquid or gas, respectively._NAis the molar transfer rate._DABis the binary diffusion coefficient. This expression derives into a formula that contemplates resistance to mass transfer due to motive force, mass transfer area and motive force or bulk due to agitation. Equation (4) stands for the volumetric mass transfer rate, whereby_NAhas units ofmol ^{m−3 s−1}.

_NA=kaΔ_CA=ka(_CAL*−_CAj)

In this study, it was necessary to consider the case where the mass transfer is in the gas–liquid interface because it represents the oxygen transfer from the bubbles to the broth. The transfer of the gaseous solute is analyzed in the same way as the transfer of liquid–liquid and solid–liquid. It is assumed that the solute is transported from the gas phase to the liquid phase. With this, Equations (5) and (6) can be deduced [3].

_NAG=_kGa(_CAG−_CAGi)

_NAL=_kLa(_CAL−_CALi)

In 1989, Chisti and Moo-Young presented a refinement of these equations by considering the following assumptions [12]: (i) Transfer through the gaseous bulk to the bubbles is fast. (ii) The resistance imposed by the liquid–gas interface is assumed to be negligible. (iii) The stagnant film around the bubbles, which is the highest transfer resistance, is minimal. This is because the tank is considered to be adequately mixed, and there is a predominance of convective diffusion. However, this film is the most considerable resistance to oxygen transfer, and the transfer in this area will be predominant. (iv) The liquid bulk is neglected because the tank is well agitated. Besides, for the case study, the viscosity is low enough not to affect oxygen transfer. This implies that the expression (6) becomes zero due to the concentration differential. Finally, Equation (7) is obtained.

_NA=_kLa(_CAL*−_CAL)

where_CALis the oxygen concentration in the broth inmol⋅^m−3and_CAL*is the oxygen concentration in equilibrium with the gaseous phase also inmol⋅^m−3known as saturation concentration or oxygen solubility in the liquid. In the case of batch cultures, the rate of oxygen consumption varies over time since the concentration of cells increases over time, and the rate of oxygen consumption is proportional to the number of cells in the culture. The biomass consumption rate, known as the specific oxygen consumption rate, also varies until a maximum is reached in the first stages of cell growth._Qois defined as the consumption rate by volume, and_qoas the specific consumption rate, depending on the biochemical cell nature and the nutritional environment. Equation (8) relates both terms.

_Qo=_qox

where x is the cell concentration; both the cell concentration and the specific consumption rate, vary for each type of culture, whether microbiological, plant cells or animal cells. Likewise, when the dissolved oxygen concentration is lower than the minimum consumption rate of the cell (critical consumption), the mass transfer depends only on the oxygen concentration in the liquid. It must be ensured that the oxygen concentration in the culture is than the critical oxygen concentration(_Ccrit)to ensure that_qois constant and independent of_CAL. Otherwise, if_CALfalls below_Ccrit,_qohas a linear tendency dependent on oxygen concentration. If it is assumed that there is no oxygen accumulation in the reactor as the cell culture consumes oxygen at a constant rate in a given period,_NAin Equation (7) can be equated to_Qoto obtain Equation (9).

_kLa(_CAL*−_CAL)=_qox

With this equation, several parameters of interest for the bioreactor design can be set up.

Minimum_kLato Cell Culture

This term refers to the minimum mass transfer coefficient(_kLa)for_CALto be higher than _Ccrit(Equation (10)).

_(kLa)Crit=_qox(_CAL*−_Ccrit)

With these parameters, it is known that the agitation in bioreactors must be slow since it needs a low shear rate to ensure that the cells do not present a risk of rupture, especially in mammalian cells. The speed ranges between 100 and 400 RPM. 2.4. Governing Equation For the equations that govern the fluid, it is essential to define three fundamental equations: continuity, momentum and energy balances. However, only the first two were considered in this study as it is assumed that the reactor is adiabatic.

2.4.1. Continuity Equation

The continuity equation (Equation (11)) models the conservation of matter in a continuous medium for any fluid.

∂_ρ ∂t+∇.(_{ρ v→} )=0

The first term describes the transfer rate of matter in the medium and the last term accounts for the mass entering and exiting the medium by diverging velocity multiplied by density.

2.4.2. Momentum Equation

The momentum equation describes the transfer of momentum through all the fluid in the system, analyzing all the possible terms that affect the fluid from the perspective of this transport phenomenon (Equation (12)).

∂_{ρ v→} ∂t+∇.(_{ρ v→ v→} )=−∇_P +∇.[μ(∇v¯+∇^vT)]+_{ρ g→} +F→

On the left-hand side, there are two terms: the first term describes the rate of change of momentum through the control system and the second considers the convection of the fluid. In contrast, on the right-hand side, the first and second terms are the molecular contributions to the phenomenon and the last two terms evaluate the possible external forces that could affect the fluid, such as gravity. The following expression can be obtained from the momentum equation, assuming constant density and viscosity (Equation (13)).

ρD_v→ Dt=−∇_P +μ^∇2 _v→ +_{ρ g→}

The definition of substantial derivative is used to accommodate the terms on the left, and external forces are neglected except for gravity. This is the Navier–Stokes (NS) equation, which is the point of departure for the flow analysis [11].

2.4.3. Turbulence Model

As the above is defined for Re lower than 10, the regime is laminar. Then, to not make modifications to the turbulence model, it was assumed that, for regimes in transition, the turbulence models described below apply. The RANS (Reynolds-averaged Navier–Stokes) model allows for solving the NS equations. Nevertheless, it requires of other models to be precise, since there are extra parameters that must be calculated depending on the system, as it is the turbulent stress which counts at the same time with the term known as Eddy’s viscosity [13].

The models for calculating this term of viscosity change depending on the number of equations that try to describe it; one of these is the so-called k-ε turbulence model [14]. The k-ε turbulence is a method of searching for the Eddy Viscosity term present in the RANS approximation, using two equations: the k that describes the kinetic energy of the system and the ε which by definition is the turbulent dissipation term and gives information on how much kinetic energy in the turbulence is converted into thermal energy [15]. The model currently used is a variation called Realizable k-ε, which is responsible for increasing the robustness of the initial model generating essential changes to understand systems with vorticity in turbulence; therefore, it is widely used in industry [15].

2.4.4. Eulerian Multiphase

To understand the Eulerian multiphase model, it is essential to understand the Eulerian framework. For multiphase flow analysis, there are two approaches, the Lagrangian and Eulerian frameworks, which do the same work modeling systems with different phases. However, both have different approximations, mostly in the way they are resolved. The first one analyzes point by point the distinct phases. In contrast, the second one approximates the phase as a continuum, and it is computed not point by point, but instead of a fixed grid where the phenomenon will develop [16].

The Eulerian model divides every dependent variable into an average (U) and a fluctuating component (u), generating the Eulerian wind field vector [u, Equation (14)]. On the other hand, ⟨ ⟩ is the average ensemble operation, thus it is correct to define c =⟨c⟩+c’, where c is the concentration of a dispersed phase and^c′ is the fluctuation part [17].

∂⟨_ci⟩∂t=−U.∇⟨_ci⟩−∇.⟨_ci’u⟩+_Dm ^∇2⟨_ci⟩+⟨_Si⟩

where D_m is the molecular diffusivity and_Siis the generation term for species, where chemical reaction must be considered, and the termsU.∇⟨_ci⟩,∇.⟨_ci’u⟩andD^∇2⟨_ci⟩ correspond to the advection, turbulent diffusion and molecular diffusion of the species [16].

This model assigns the momentum and continuity equations in phases and then joins them by relations of pressure and the interphase exchange coefficients, according to the phase-type. This model is widely used for bubble and fluidized bed column systems, which is why it was used to perform the aeration in the bioreactor [16].

In addition, the Volume of Fluid model (VOF) works fine for systems with large structure phases with a small total contact area. This model is used because of its numerical efficiency, even though its risky approximation due to the behavior of the gas. However, since the system is not modeled to analyze the bubbles, instead, it determines the functioning of the multiphase flow, considered a reasonable estimate [18].

3. Materials and Methods 3.1. Experimental Methods

3.1.1. Impeller Design

In this study, a characterization of a bench-top bioreactor was carried out based on power number, flow patterns and oxygen transfer rate. Three types of impellers were characterized: Rushton, 3-paddles, and propeller (Figure 1A). They were designed according to the dimensions of a standard tank, which have in general a diameter of 0.06 m [19]. Additionally, to make comparisons in the performance of the diameter of the impellers, a new impeller with a smaller diameter was designed. Lastly, different impeller combinations were studied, as shown in Figure 1B. The CFD software STAR-CCM+ v13.04 (SIEMENS) was used to simulate parameters such as power, torque and flow patterns of the bioreactor. The impellers’ dimensions were recreated in the Autodesk Inventor Software® v.2018. The choice of this software was due to the licenses available when this research project was carried out. Any 3D modeling program (3D CAD) can be used since it can export files in Parasolid format (x_b or x_t). Besides, the modeling assistant incorporated in the STAR-CCM+ software can be used for this purpose. The isometric planes of each impeller and the vessel of the bioreactor are presented in Figure A1.

3.1.2. Power Curve Determination

It was necessary to determine the Reynolds values and power numbers (_Np) for each flow regime to create the power curve for each impeller. This is essential to establish the behavior of each impeller from the laminar to the turbulent regime. The Reynolds number, as shown in Equations (1) and (2), depends onρ, N, D and μ. The power number, on the other hand, depends onP, ρ, N and D. The power curve was calculated, ranging from10<Re<100,000.

Additionally, some variables were fixed, such as agitation velocity and impeller diameter, as shown in Table 1. The experiments were performed at a constant temperature of 15 °C (Bogotá’s ambient temperature). However, as previously mentioned, it was necessary to adjust the fluid’s viscosity to precisely determine the power curve. For this, a viscosity curve of a glycerin–water solution at different concentrations (%v/v) was carried out, as shown in Figure A2. Equation (15) is obtained to describe the change of the broth viscosity as a function of the water percentage in a water–glycerin solution.

μ=2849.1⋅_xw−3.138

3.1.3._kLaDetermination by Gassing Out Method

The oxygen content in a tank, at any given time, is a result of the oxygen consumption by the cells within the culture medium and the continuous oxygen transfer to the medium. The oxygen transfer rate to the tank must be considered to determine the oxygen uptake rate (OUR) of the cells. For this, it was necessary to calculate the oxygen transfer rate coefficient(_kLa), by controlling parameters such as the airflow (in VVM, the volume of air per volume of medium per minute). Thus, oxygen at the system outlet is considered to be the oxygen concentration fed into the system minus the consumption of the cells, assuming that the dissolved oxygen is only used by the cells according to their OUR. Once the bioreactor has been assembled, the pH electrode and the DO probe were calibrated. The latter was set at a concentration of 20% initially. Subsequently, the DO was varied according to an established experimental design (described below).

Moreover, the effects of the impellers’ combinations (Figure 1B) in the_kLawere considered. The oxygen in the reactor was displaced with nitrogen using the gassing out method. This consisted of nitrogen or carbon dioxide injection until a minimum critical DO percentage was reached (for this study, nitrogen was used). After this limit was achieved, the air valve was opened to reach the required value.

The concentration of DO was measured every 15 s using an oxygen probe until it stabilized. Finally, by using Equation (16), which is a treatment of Equation (9), the term q_ox can be replaced by a derivative definition (dC_AL/dt). At the same time, the differential term can be integrated between two points (t₁ and t₂), resulting in Equation (16). The slope of linear regression was determined, corresponding to the_kLa.

(_t2−_t1)_kLa=ln(1−_CALCAL*)

where_ticorresponds with the time and has units ofs. The value for oxygen solubility(_CAL*) was obtained using Equation (17), initially proposed by Doran [3].

_CAL*=14.616−0.3943 T+0.007717 ^T2−0000646 ^T3

where_CAL*has units ofg ^l−1and temperature (T) has units of° C. It is necessary to clarify that a common problem in this type of analysis is the response time of the DO probe. However, taking into account the nature of the controller (P&ID), the manufacturers adjusted the proportional constants so that the sensor response was immediate [20]. Additionally, studies such as the one carried out by Spanjers and Olsson show that making a first-order adjustment to the curve has good results because it takes the area in which the value of the slope changes in a linear way and thus the response time does not significantly affect the results obtained, and this is one of the assumptions made when finding the_kLa values in this study [21].

3.1.4. Experimental Design to Establish Operational Effects in Dissolved Oxygen

Since it is desirable to determine the operating conditions that affect OTR, an experimental design with three main factors that can influence bubble retention and mass transfer within the bioreactor was proposed. Table 2 shows the factors and levels evaluated in this study.

Once all factors and levels were established, the statistic software Minitab^® (version 18, Minitab LLC, State College, PA, USA) was used to create the experimental design order shown in Table A1 and to statistically evaluate the data.

3.1.5. Flow Patterns

The importance of the flow patterns resides in the capacity to maintain the oxygen bubbles as long as possible, avoiding coalescing, which decreases the rate of oxygen transfer to the culture medium in stirred tanks. In addition to oxygen transfer and power calculation, flow patterns depend on aspects such as the agitation velocity, the fluid properties such as viscosity and the impeller type. Simulations were carried out to define the flow patterns with each impeller. 3.2. CAD and Mesh Construction

3.2.1. Bioreactor Dimension and Geometry Design in Autodesk Inventor Software

To carry out the corresponding simulations and to be able to confirm with the experimental data, it was necessary to determine the dimensions of the reactor studied. The measurements were taken for the geometry of the Brunswick Bioflo/CelliGen 115 bioreactor (Eppendorf, Enfield, CT, USA), equipped with temperature, pH and DO probes. The pH is sensed by a gel-filled pH probe in the range of 2.00–14.00. Control is maintained by a P&I controller which operates peristaltic pumps, assigned to perform acid or base addition, or which controls the use of gas (es) for this purpose. A Polarographic DO electrode senses the dissolved oxygen, and a P&I controller maintains control by changing the speed of agitation, controlled in the range of 0–200%. It is thermal mass flow controller-regulated flow rate, and the percentage of oxygen in aeration with a digital display in 0.1% increments [20], stirring system and anti-foam control, with a capacity of 5 L [22]. The final model is shown in Figure 1C.

3.2.2. Mesh Independence and Preliminary Configuration

It was proposed to determine mesh independence, both quantitative and qualitative, to obtain the best results for the simulations. For both techniques, a base value was chosen in the mesh. This value was determined based on the conditions of the geometry. The recommended amount is a deviation of 30%, as defined by Celik et al. [23]. The next thing was to choose the variable on which the comparison between meshes would be made. The power and shear rate were determined since they are the primary response variables analyzed in this case study. Simulations were performed, and the results are interpreted in two ways. The first is a qualitative technique in which visual aids were used. In contrast, the second is only quantitative in which the Grid Convergence Index (GCI) was determined, whose process is described below [23].

A base mesh was made for the bioreactor. This mesh configuration was called “normal” and had four prism layers to resolve the boundary layer near the wall correctly. Table 3 presents the mesh configuration used in this study.

Once the base mesh was defined, the next step was to create the different meshes. The first was called “coarse” since it had a smaller number of cells, while the second was called “fine” because it had a larger number of cells. Each simulation was performed using the same physical parameters.

Qualitative Method

This technique was required to graph the chosen variables, in this case, power and shear rate, and to analyze the trend they have concerning mesh size. The ideal situation is to obtain a constant trend as the base size decreases. However, this technique does not consider the extrapolated or interpolated values. Therefore, many times, in between meshes must achieve the expected asymptotic result [24]. The results obtained in the independence were shown through the different reports generated.

Quantitative Method–GCI Calculation

The GCI was used to analyze the convergence factor that a mesh has for distinct types of variations, as, in this case, for the base size. The step-by-step approach to calculating this index and correctly performing mesh independence, is described in Appendix A, according to the fluid engineering division of the American Association of Mechanical Engineers (ASME) [24]. With the GCI, it is possible to determine whether a subsequent refinement of 30% less was necessary, giving a percentage that must be carefully examined since it is up to the researcher’s consensus whether to continue with new simulations or not. For this case, an acceptability percentage of 10% was defined, since each of the fine meshes generated has many cells, so that the time of convergence would indiscriminately increase if there were a lower criterion.

3.3. Modeling Approach The modeling was focused on the power and flow patterns analysis of the bioreactor. Therefore, two different routes, with distinct initial and boundary conditions, were considered. The simulations were performed on two virtual machines (Microsoft, Redmond, WA, USA) with 12 cores and 120 GB of RAM provided by the Information and Technology Services Management (Dirección de Servicios de Información y Tecnología, DSIT) of the Universidad de Los Andes.

3.3.1. Power Analysis

The power analysis was performed as a single-phase steady-state simulation using water as a fluid that fills the bioreactor, varying the revolutions per minute to modify the flow regime of the system. The only condition of the process was the moving reference frame (MRF), as shown in Figure 2. The convergence criteria used for the power analysis simulation was to inspect the residuals for the different solvers until they reached a value of less than 1 × 10⁻⁴ when the variations between the response were minimal.

3.3.2. Flow Patterns Analysis

On the other hand, an implicit unsteady gas–liquid simulation was conducted in other to analyze the flow patterns in the bioreactor. A Rigid Body Motion (RBM) scheme was implemented to predict the flow patterns with a rotation involved. The air is introduced into the system through a perforated probe, leaving an opening located in the cap of the bioreactor (Figure 2). Using the VOF scheme, the system initializes by filling the tank with water up to 3/4 of its height, and the rest is completed with air. By simulating in the transient state was essential to determine the duration of the physical time, which, in this case, was aimed to analyze the system until the impeller made three revolutions. The flow becomes quasi-static from the second 0.0002, where the power of the tank stabilizes and remains constant. It takes 0.72 s for the agitator to make the three expected rotations.

4. Results and Discussion 4.1. Mesh Independence

Due to the nature of this study, it was necessary to determine the quality of the simulated models, so a mesh independence analysis was needed. Table 4 shows the different errors and GCI. The absolute error only considers the interior interval between the meshes already obtained, so it is an indicator of how much the results vary between the fine mesh and the normal mesh. The extrapolated error gives information on how precise the mesh is for base size values smaller than fine; thus, it provides information beyond that which could be simulated.

For Brunswick Bioflo/CelliGen 115 bioreactor, high values of both GCI and errors were observed. These values imply that this model has an oscillatory convergence. Therefore, refinement was necessary until it reached the desired asymptotic trend. However, this method did not consider computational time. Hence, its refinement was carried out until the computational expense and time allowed it was met.

Mesh Independence Analysis

The GCI values were analyzed, considering an acceptability value of 10% in computational time. The meshes defined across the qualitative method were used, bearing in mind that, even if the quantitative approach is precise, it does not consider the computational time required to carry out the simulations. Figure 3 shows the power and shear rate versus the number of cells. It can be observed that, as finer meshes (higher numbers of cells) are used, the evaluated variables tend to keep their values unchanged, which is expected behavior. It should also be noted that, as the number of mesh cells increases, the computational time also increases. As for the shear rate, the behavior between the different meshes was not constant, but the variation from normal to fine mesh was minimal.

It was necessary to add a new mesh between normal and fine, which was called semi-fine (reducing the base size by 15%), to validate if there was a consistent behavior in that area. The results of mesh independence for the semi-fine mesh (selected mesh) can be seen in Table 5. It was considered that the computational time for a finer mesh would generate an increase in the cost of each simulation, thus the semi-fine mesh was chosen as the main mesh for the different simulations.

4.2. Simulation Model Validation Initially, the analysis was realized with only water, evaluating the torque for different angular velocities in the Brunswick Bioflo/CelliGen 115 bioreactor, using various impellers. However, it was not possible to measure the torque for low RPM due to sensor sensitivity. Nevertheless, it was possible to obtain a single point at 100 RPM for inclined paddles. Moreover, it was necessary to design a script in JAVA^®, which only uses the geometry and mesh of the bioreactor as inputs, to calculate the Re vs. Np curve automatically.

The inclined paddles were the first approach used to validate the simulated model with the experimentation, as shown in Figure 4A. Notice that the experimental curve fits well at the beginning of the simulated power. However, this tendency changes at higher RPM, where the experimental data end at 1.5 W, and the simulated ones increase to approximately 1.9 W. This variation between the experimental results and the CFD model could be due mainly to a couple of factors: the first one is purely experimental, because the torque sensor worked better for more viscous substances. The second one is due to the CFD model. Therefore, we propose in future studies switching the RANS model, to observe if there is any significant variation in results (Figure 4A).

On the other hand, it is possible to create a new plot analyzing the behavior between the experimental and simulated power, linking the data with RPM in Figure 4B. As shown in Figure 4B, it is reasonable to say that both power values increase as RPM increase, and the simulation model behaves correctly at low RPM. However, it is impossible to determine the accuracy of the experimental model without the analysis at low RPM. As mentioned above, the solution for this issue was to create a new experimental model where the RPM variable is changed with the Re number, using different viscosities. The purpose of this is to be able to use the definition of Re (1) and change the viscosity according to Equation (15). This facilitates the construction of theRe vs. _Np curve, as shown in Figure 5. This is how Figure 5A was obtained, which shows a similar linear behavior of inclined paddles for laminar and transition regime for the experimental and simulated curves. A continuous line for turbulent regions is also observed.

TheRe vs. _Np curve for Rushton is shown in Figure 5B. Five experimental points were used for comparison. Meanwhile, for the simulated curve, additional points were obtained. As expected, both curves have a similar trend. However, the last experimental point for_Npcan be considered an experimental error. Nevertheless, the differences between both curves are minimal or negligible with relative errors of about 28%, which is still consistent due to possible experimental errors affecting the measurements.

A small propeller was designed, 3D printed, and experimentally tested to validate the model. If the results show at least similar trends, then the simulation can be regarded as adequate. The results are shown in Figure 5C. The behavior and performance of the experimental data displayed that the model presents some small variations, usually of experimental errors. However, due to the closeness of the data, it can be asserted with certainty that the model is accurate.

4.3. Power Analysis

Once the simulation model was confirmed, the power curves for all the impellers were reproduced using CFD models. It is necessary to clarify that, experimentally, the torque of the inclined paddle, Rushton, and small propeller agitators was measured. The scope of the model extends to cases where there is more than one agitator in the reactor, thus theRevs._Np curve of the proposed combinations (Figure 1B) was also determined via simulation.

As can be seen in Figure 5A, the experimental data and the simulation data do not have the same range due to measurement difficulties atRelower than 700. One can argue that the inclined paddles have a low-power consumption behavior for the laminar regime. Contrarily, the tendency for turbulent and transition regimes is of high-power consumption.

The generic propeller of the bioreactor was only simulated, as shown in Figure 5D. With the CFD model approach, more data can be collected, even for regimes that cannot be replicated in the laboratory. The generic propeller presents a high-power consumption at the laminar and transition regime. On the other hand, the small designed propeller has a lower power consumption. For Re equal to 10, the propeller has an Np value of 323.45, while the small designed propeller has an Np value of 65.30, a difference of an order of magnitude. Furthermore, for turbulent flows, the propeller has a better performance than the designed impeller.

As previously demonstrated, the model presented an excellent performance. Different combinations of impellers were analyzed, taking advantage of the model, determining and predicting improvements in the stirring process. Figure 5E shows, for the laminar regime, that the differences between propeller–Rushton and Rushton–paddles are negligible. Similarly, the differences between the last two and the paddles–propeller are around one order of magnitude, having similar behaviors. On the transition regime, the data diverge since the paddles–propeller has an order of magnitude lower than the other two combinations. In contrast, in the turbulent regime, the responses of Rushton–paddles and paddles–propeller are similar, whereas the propeller–Rushton has a different performance, with a lower order of magnitude.

4.4. Flow Patterns

Following the bioreactor distinction, it was decided to model in CFD the behavior of the flow patterns of the primary impellers (Figure 6). It was possible to understand the importance of agitation to increase the bubble residence time in the medium, promoting mass transfer.

In Figure 6A, the axial flow in the tank is visible, where the fluid is propelled parallel to the impeller axis [26]. This is important since these patterns have an upward tendency, crashing and returning by the action of the hydrodynamic force, causing vortices to form at the tips of the impeller. Then, it is possible to notice that the fluid has a remarkable perturbation over the whole tank, which is favorable because the bubbles found in the medium can effectively reach the cell culture, improving diffusion. Likewise, due to the size and shape of the blades, oxygen may remain longer in the system.

On the other hand, the Rushton impeller has a radial behavior due to the shape and geometry of its blades. In Figure 6B, the impeller moves the fluid to the tank walls, where it collides and is propelled to the bottom and the top of the tank, a typical behavior in radial flows, as exposed by Spogis [26]. It is important to note that, due to its nature and size, a Rushton impeller is not adequate at keeping air bubbles in the medium because it concentrates only on one specific area of the entire bioreactor. This phenomenon makes the Rushton–Rushton configuration the one typically used to prevent the formation of dead zones.

The propeller has an axial behavior, which is similar to the inclined paddles. In this case, the flow goes in the rotational direction, causing the fluid to hit the bottom of the tank and return, forming a vortex of greater magnitude than from inclined paddles around the blades. Even if this impeller is smaller than the paddles, it can impart a substantial momentum over the entire fluid column, something that Rushton does not accomplish. Agitation is vital in oxygen diffusion. Depending on where the sensor is located in the vessel, there may be errors in the measurement. Nevertheless, it is assumed that perfect agitation is achieved and that the oxygen diffusion is almost instantaneous along with the tank. The most crucial role of the agitators and, in particular the flow patterns formed, is how fast the oxygen diffuses into the medium. On the other hand, the flow patterns strongly define the movement of the bubbles; with this, it is understood that radial patterns promote the time in which the bubbles remain in the reactor since a type of barrier can be formed that prevents the bubbles from moving in a vertical direction towards the air outlet. 4.5. Oxygen Diffusion

4.5.1. Experimental Design

Considering the factors and levels proposed in the experimental design (Table 2), the measurement of each experimental configuration was carried out evaluating the concentration of dissolved oxygen in the bioreactor. The statistical results are shown in Table A1. Likewise, an Analysis of Variance (ANOVA) was carried out to determine each variable statistical significance and the interactions between those variables Table A2 shows the results for a significance of 5% (0.05). As shown in Figure A3, the average of the response variable increases as the airflow and stirring speed increase. This is expected since, at more significant airflow, more bubbles will be present in the medium, promoting mass transfer. Concurrently, increasing the stirring velocity stimulates the turbulence of the system, increasing the retention of bubbles in the medium. The type of impeller also has a significant effect on the response variable. However, the lowest average is obtained when a propeller-type impeller was used. The difference between the paddle impeller and the Rushton impeller is approximately 5%, being higher than the average obtained when using a Rushton type impeller. According to the ANOVA, each of the variables analyzed is statistically significant since the type of impeller, the airflow and the agitation velocity have p-value values of 0.002, 0.004 and 0.004, respectively. It is important to note that the change in agitation velocity has a more substantial impact on the average of the response variable. If each factor were evaluated separately, the maximum possible average would be achieved by setting the stirring speed to 250 RPM. Considering only the graph of main effects (Figure A3), the configuration that ensures a higher rate of transfer of oxygen is a Rushton-type impeller, an airflow of 5 L/min and 250 RPM.

As for the interactions, only binary interactions were studied. The results obtained show that the only significant interaction (p-value = 0.034) is the one between the type of impeller and the agitation velocity. This indicates, as already mentioned, that these two factors directly affect the mean DO and are following the results obtained in Figure A3. The interaction between the type of impeller and any other factor influences the average response differently than the main factors Figure A3. As the impeller speed increases, the performance of the impeller is affected. However, once the stirring speed is increased, the paddle-type impeller presents better results in DO. When blades and a stirring speed of 250 RPM were used, the highest possible value was achieved in the response variable. With the propeller-type impeller, it is observed that when the agitation velocity is increased, both propellers have the same performance. With the interaction between the airflow and the impeller type, the same pattern is obtained, being better the performance of the paddles when there is more significant airflow. Finally, concerning the interaction of the air with the agitation velocity, this is consistent with Figure A3.

4.5.2._kLaDetermination

Once all the bioreactor–impeller configurations (single or combined) were obtained, the mass transfer coefficients(_kLa)were determined by the gassing out method. Initially, the measurements with one impeller were carried out, followed by the impeller combinations measurements, at four cases of operation for airflow and agitation velocity: (i)2.5 l ^min−1 and 100 RPM; (ii)2.5 l ^min−1 and 200 RPM; (iii)5 l ^min−1 and 100 RPM; and (iv)5 l ^min−1 and 250 RPM.

One Impeller

For the case where each impeller was evaluated separately, measurements were made for the four proposed impellers (Figure A4). Additionally, in Table 6, the exact slope values of each curve are given, corresponding to the_kLafor each configuration, together with its corresponding power and shear rate values.

Different analyses were obtained, as shown in Figure A4. First, the effect of agitation on propeller performance: as agitation increases, the slope of the small propeller was higher than the standard propeller. For cases where great revolutions are used, it is preferable to use the small propeller because it requires less power, and the shear rate is lower. Therefore, there is less risk of cell wall rupture. This reflects, among other things, the effect of diameter on power consumption and shear rate, by comparing the runs made with the propeller (Runs 5–8) and the small propeller (Runs 9–12). Specifically, when comparing Runs 5 and 9 where the conditions are the same, the_kLavalue obtained with the small propeller is lower by only 0.0005^s−1, but, when observing the power consumption, it is lower by one order of magnitude. Likewise, the shear rate is lower for the case of the small propeller (2.376 against 0.723^s−1). This indicates that the reduction of the diameter is a factor to take into account since the same results can be achieved with a standard propeller of 0.06 m with lower power and less shear stress on the cells. This analysis can be performed by comparing any of the runs between the large propeller and the smaller diameter propeller.

Depending on the case, the Rushton impeller is better than the paddle. In most cases, both impellers present favorable results for DO. However, at low velocities and low airflow, the difference is more remarkable, and Rushton has the best performance. In all other cases, the choice of the impeller will depend on the power required and the shear rate of each one.

Impeller Combinations

Special consideration was made for the impeller combinations, that is the location of each one. The possible combinations are shown in Figure A5, and it was decided to measure the concentration of dissolved oxygen for cases where the positions of the impellers are inverted. This is only for situations where there is an airflow rate of 5l ^min−1, because, as previously shown, the variation in the speed of agitation has a more significant effect on the diffusion of oxygen. Therefore, it was assumed that inverting the positions does not affect the power or the shear rate.

As shown in Figure A5, in general, the combination of Rushton–paddles has the most significant slope in most cases. However, when the stirring speed increases to 250 RPM, the inverted propeller–paddle combination has the best performance; in turn, this combination at operating conditions of 5l ^min−1and 250 RPM has the highest slope of all, which means they have the highest overall_kLavalue. This result is remarkable because the performance of the propeller is lower than the rest, even though the combination of Rushton–propeller impellers had the lowest slope in all cases. It is important to note that flow patterns play a vital role in these measurements. With the mixing of two impellers, the flow patterns are also combined. Placing the propeller at the bottom of the tank cuts dead zones at the bottom of the bioreactor, and the blades increase turbulence. This achieves the highest efficiency in transferring oxygen from the bubbles into the broth. Finally, to give applicability to the_kLa values obtained for each impeller and combination to different operating conditions, these values are presented in Table 6. For this study, only one nozzle size was used for_kLa determination. However, this is an important factor when optimizing the oxygen transfer rate to the culture. Studies show that the smaller is the nozzle size, the smaller are the bubbles, and the higher is the oxygen transfer rate. This is since the smaller is the bubble, the longer is the bubble retention time in the bioreactor. Likewise, as the volume area ratio increases, the smaller is the bubble size, the larger is the surface area and the higher is the oxygen diffusion rate [27,28].

Shear Rate

To analyze the data obtained from the shear rate for each of the operating conditions proposed in Table 6, it was decided to compare these with maximum shear values of common cell strains used in biotechnology. Data were taken from seven different cell strains: the bacterium Escherichia coli (EC), the yeast Saccharomyces cerevisiae (SC) [29], the mammalian Chinese hamster ovary (CHO) cell [30,31], Aspergillus glaucus (AG) [32], Trichoplusia ni (NT-368) and Spodoptera frugiperda (SF-9) [33] and HeLa cells [34,35]. The properties of the cells and their cultures, such as shear rate and viscosity, are presented in Table A3. With this information and the shear rate values in Table 6, Figure 7 was obtained, which compares each of the operating conditions studied with the typical shear values of the cells, to determine which is the most appropriate for a selected range of cultures.

In Figure 7, it is observed that the shear rate values obtained for each operating condition are below the maximum shear rate values of the EC and SC cultures. Likewise, all values exceed the maximum shear rate value for CHO cultures. In general, cell lines such as CHO cells (mammalian cells) do not require mechanical agitation, and, in cases where agitation is present, it is usually gentle [30]. According to Figure 7, for cells such as AG, some of the conditions studied show a higher shear rate than the maximum supported by this cell type. Examples of this are Runs 2, 4 and 19. However, 26 of the 34 conditions are feasible, and it is possible to find the optimal configuration with which to grow an AG culture. For the other cell types analyzed, the maximum rate is higher than the topmost shear rate data presented in Table 6. Therefore, any evaluated operating condition is operable for cells such as EC, SC, TN-368, SF-9 and HeLa. The choice of the best agitator will then depend on the_kLavalue required to keep the culture conditions and the power needed for the agitation process, and it is ensured with the 34 operating conditions studied.

5. Conclusions

The three factors studied, impeller type, agitation rate and airflow, have a direct influence on the oxygen transfer coefficient_kLa. The agitation velocity proved to have a more significant impact on the mass transfer because, as it increases, the diffusion of oxygen in the bioreactor becomes more favorable. This is due to the turbulence increase, and the bubbles’ residence time inside the fluid is longer. The type of impeller with the best performance is the Rushton type. However, with increasing agitation velocity and the addition of constant airflow, inclined paddles are highlighted if only one impeller is used in the reactor. Another aspect that is relevant to the study, even though it was not profoundly analyzed, is the size of the bubbles. When the bubble size is smaller, the surface area–volume ratio is larger, and therefore there is a higher oxygen transfer rate. Similarly, when the bubble size is reduced, it tends to remain in the vessel for more extended periods.

The addition of a second impeller presents benefits for the oxygen transfer in mixing processes. When the two impellers are of a different type, their flow patterns are combined to provide a longer bubble residence time. Just as in the case where there is only one impeller, the variables of agitation velocity and airflow can increase the_kLavalue. The performance of the propeller-type impeller rises dramatically as the agitation rate increases. Thus, when it is decided to operate the reactor with a stirring speed of 100 RPM, the combination of paddle-type impellers and Rushton ensures more significant mass transfer in less time. However, when the reactor is operated at 250 RPM, the propeller paddle combination performs better. Concretely, the location of the propeller at the bottom and the blades at the top of the reactor so that the flow patterns favor the transfer.

Moreover, considering only the power requirement for the different impellers, it can be concluded that the propeller requires less energy to break the inertia for turbulent and transitional flows. In contrast, for laminar flows, the inclined paddles require less power. However, analyzing the combined impellers, it can be appreciated that all have an unfortunate development for laminar zones, while, for transition and turbulence, present results similar to the propeller. The propeller and inclined paddles showed flow patterns that can increase the residence time of the bubble plus add vorticity to the entire system, which makes the bulk reduced; this is important because these two impellers have the lowest energy requirement compared to Rushton. The studied conditions proved promising for bacterial and yeast cultures, but the shear rates showed a deleterious effect on mammalian cells (CHO cells). This study is of interest to the bioprocess engineering area, especially in bioreactor design, since it shows a fast, reliable and affordable alternative to ensure proper mixing and oxygen availability in the system, parameters directly affecting cell growth, yield and productivity in a bioprocess.

Follow up studies will consider improving the power analysis, changing to a multiphase system, first using the VOF equations and then moving to a purely Eulerian scheme, where the drag forces of the system are considered, analyzing whether the power is affected by the aeration of the system to reduce the model error. Each phase can be differentiated with a separate turbulence equation, handling the continuous phase with a k-omega model and the dispersed phase with a k-epsilon model, as proposed by Karpinska and Bridgeman [36,37], who presented excellent results in agitated systems with air injection. On the other hand, the present study considered the case where the temperature does not change over time. This is an assumption that may not be correct in all cases. Therefore, studying the influence of heat and mass transfer may result in a complete model from the computational perspective. This makes the model more robust, requiring additional computational resources due to the larger number of equations needed.

Figure 1.(A) Impellers used for this study (from left to right): Rushton (six-blade turbine), three-paddle helix impeller, Brunswick Bioflo/CelliGen 115 bioreactor default propeller and small propeller (abbreviated as S. Propeller). (B) The impeller combinations used (from left to right): Rushton-paddles, propellers-paddles, and Rushton-propeller. (C) The 3/4 section view 3D CAD models of: Brunswick Bioflo/CelliGen 115 bioreactor (left); and transparent view (right).

Figure 2.(A) Boundary conditions for the power analysis using a Moving Reference Frame (MRF) scheme. (B) Boundary conditions for the flow patterns analysis, with the inlet and outlet of air and the Rigid Body Motion (RBM) scheme. In both cases, the yellow area corresponds to the boundary location of the rotation zone.

Figure 3. Mesh independence for: (A) power (●) and shear rate (■); and (B) computational time on Brunswick Bioflo/CelliGen 115 bioreactor.

Figure 4.(A) The first experimental process approximation for the calculation of power on Brunswick Bioflo/CelliGen 115 bioreactor using the impellers with inclined paddles. The curves for simulated power (●) and experimental data (■). (B) Analysis between the experimental and simulated data for the power of Brunswick Bioflo/CelliGen 115 bioreactor using impellers with inclined paddles.

Figure 5. Power number curve for the different regimes in Brunswick Bioflo/CelliGen 115 bioreactor using inclined: paddles impeller (A); Rushton impeller (B); small designed propeller (C); and default propeller (D). Curve for simulated power number, experimental data, and literature data [25,26]. (E) The simulated power number for different Re numbers using various impeller combinations can be observed: propeller-Rushton, Rushton-paddles, and paddles-propeller.

Figure 6. Flow patterns for Brunswick Bioflo/CelliGen 115 bioreactor at 250 RPM and 2.5l min-1for: (A) inclined paddles; (B) Rushton; and (C) propeller.

Figure 7. Shear rate values for each operation condition compared with the maximum shear rate of three different cell cultures.

Impeller Type	Diameter (cm)
Rushton	6.09
Paddles	11.46
Propeller	10.40
S. Propeller	6.01

Factors	Levels
Impeller type	Propeller
	Rushton
	Paddles
	Small Propeller
Airflow (L/min)	2.5
	5
Agitation velocity (RPM)	100
	200

Parameter	Brunswick Bioflo/CelliGen 115 Bioreactor
Base Size (cm)	2.55
Relative target size (to base size) (%)	10
Relative minimum size (to base size) (%)	10
Relative prism layer total thickens (to base size) (%)	10
Number of prism layers	4

GCI Parameters		Value
ϕ= Power	ea21 (%)	3.98
	eext21(%)	4.85
	GCI21 (%)	6.38
ϕ= Shear rate	ea21 (%)	0.66
	eext21(%)	0.27
	GCI21 (%)	0.34

Parameter	Brunswick Bioflo/CelliGen 115 Bioreactor
Mesh	Semi-fine
Number of Cells	2.58 × 106
Angular Velocity (RPM)	600
Power (W)	44.59
Shear Rate (s−1)	83.80

Impeller	(l min−1)	(RPM)	Run ID	kLa(s−1)	kLa(h−1)	Power (W)	Shear Rate (s−1)
Paddle	5.0	100	1	0.0025	9	0.192	15.211
Paddle	5.0	250	2	0.0061	21.96	7.714	32.857
Paddle	2.5	100	3	0.0017	6.12	0.192	15.211
Paddle	2.5	250	4	0.0045	16.2	7.714	32.857
Propeller	5.0	100	5	0.0026	9.36	0.069	2.376
Propeller	5.0	250	6	0.0029	10.44	0.940	8.631
Propeller	2.5	100	7	0.0019	6.84	0.069	2.376
Propeller	2.5	250	8	0.0018	6.48	0.940	8.631
Small Propeller	5.0	100	9	0.0021	7.56	0.001	0.723
Small Propeller	5.0	250	10	0.0032	11.52	0.013	3.948
Small Propeller	2.5	100	11	0.0013	4.68	0.001	0.723
Small Propeller	2.5	250	12	0.0022	7.92	0.013	3.948
Rushton	5.0	100	13	0.0032	11.52	0.020	1.306
Rushton	5.0	250	14	0.0056	20.16	0.334	8.446
Rushton	2.5	100	15	0.0023	8.28	0.020	1.306
Rushton	2.5	250	16	0.0044	15.84	0.334	8.446
Paddle–Propeller	5.0	100	17	0.0032	11.52	0.178	13.659
Paddle–Propeller (inv)	5.0	100	18	0.0031	11.16	0.178	13.659
Paddle–Propeller	5.0	250	19	0.0062	22.32	2.869	37.119
Paddle–Propeller (inv)	5.0	250	20	0.0080	28.8	2.869	37.119
Paddle–Propeller	2.5	100	21	0.0019	6.84	0.178	13.659
Paddle–Propeller	2.5	250	22	0.0049	17.64	2.869	37.119
Paddle–Rushton	5.0	100	23	0.0034	12.24	0.231	14.968
Paddle–Rushton (inv)	5.0	100	24	0.0031	11.16	0.231	14.968
Paddle–Rushton	5.0	250	25	0.0056	20.16	3.702	41.458
Paddle–Rushton (inv)	5.0	250	26	0.0047	16.92	3.702	41.458
Paddle–Rushton	2.5	100	27	0.0020	7.2	0.231	14.968
Paddle–Rushton	2.5	250	28	0.0056	20.16	3.702	41.458
Propeller–Rushton	5.0	100	29	0.0030	10.8	0.069	3.727
Propeller–Rushton (inv)	5.0	100	30	0.0026	9.36	0.069	3.727
Propeller–Rushton	5.0	250	31	0.0028	10.08	1.111	13.792
Propeller–Rushton (inv)	5.0	250	32	0.0043	15.48	1.111	13.792
Propeller–Rushton	2.5	100	33	0.0018	6.48	0.069	3.727
Propeller–Rushton	2.5	250	34	0.0028	10.08	1.111	13.792

Author Contributions

Conceptualization, C.G.D., D.A.C.L., N.R. and L.H.R.; methodology, N.R. and L.H.R.; software, L.A.R., E.L.P. and N.R.; validation, L.A.R., E.L.P. and L.H.R.; formal analysis, L.A.R., E.L.P., N.R. and L.H.R.; investigation, L.A.R. and E.L.P.; resources, N.R. and L.H.R.; writing-original draft preparation, L.A.R. and E.L.P.; writing-review and editing, L.A.R., E.L.P., C.G.D., D.A.C.L., N.R. and L.H.R.; visualization, L.A.R. and E.L.P.; supervision, N.R. and L.H.R.; project administration, N.R. and L.H.R.; and funding acquisition, N.R. and L.H.R. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Universidad de Los Andes-Fondo de Apoyo para Profesores Asistentes (L.H.R, FAPA PR.3.2017.4627).

Acknowledgments

The authors would like to express their gratitude towards the information and technology services management (Dirección de Servicios de Información y Tecnología, DSIT) of Universidad de Los Andes for facilitating the computational resources for the project and their technical support.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

Romanic symbols

a	[m2]	Gas-liquid interfacial area
CAG	[mmoll-1]	Oxygen concentration in the gaseous phase
CAL	[mmoll-1]	Oxygen concentration in the liquid phase
CAL*	[mmoll-1]	Oxygen concentration in equilibrium with gaseous phase (oxygen solubility)
Ccrit	[mmoll-1]	Critical oxygen concentration to ensure the cell culture growth
D	[m]	Diameter
e21	[-]	GCI relative error-index
DF	[-]	Degrees of freedom
DAB	[m2 s-1]	Binary diffusion coefficient
Dm	[m2 s-1]	Molecular diffusivity
GCI	[-]	Grid Convergence Index
JA	[mols-1m-2]	Molar flux of componentA
kG	[s-1]	Mass transfer coefficient in the gaseous phase
ki	[s-1]	Mass transfer coefficient ofi
kL	[s-1]	Mass transfer coefficient in the liquid phase
MS	[-]	Mean square
N	[RPS]	Angular velocity
NA	[mols-1]	Molar transfer rate of A
Np	[-]	Power number
Nq	[-]	Pumping number
P	[W]	Power
qo	[molg-1s-1]	Specific uptake rate
Q	[m3 s-1]	Total flow
Qo	[moll-1s-1]	Volumetric oxygen uptake rate
Re	[-]	Reynolds number
RANS	[-]	Reynolds-Average Navier-Stokes
SS	[-]	Sum of square
t	[s]	Time
T	[°C]	Temperature
v→	[m s-1]	Velocity vector
VOF	[-]	Volume of Fraction
x	[gl-1]	Cell concentration in the broth
xw	[%]	Water fraction in solution
y	[-]	Cartesian plane coordinate

Greek symbols

μ	[Pas]	Viscosity
ρ	[kgm-3]	Density
ϕ	[-]	GCI analysis variable

Appendix A. Data S1. GCI Step by Step Calculation

Step 1. A representative size is defined for the mesh; this value can be the base size, minimum size or target size, among others. In the investigation, it was defined that the most representative value for the mesh is the base size since all other parameters can depend on it.

Step 2. Three series of meshes are selected to perform the simulations in such a way that the value defined in Step 1 has a variation of 30% above and below its value. In addition, it is chosen which critical variable (ϕ) is going to be evaluated, in this case, the power and the shear rate. With two different variables, it is necessary to carry out two different analyses.

Step 3. With h1 < h2 < h3 and r21 = h2/h1, r32 = h3/h2, the apparent order p of the method is calculated using the following set of equations:

p=1ln(_r21)|ln|_ε32/_ε21|+q(p)|

q(p)=ln(_r'21p-s_r'32p-s)

s=1·sgn(_ε32/_ε21)

where

_ε32=_ϕ3-_ϕ2

_ε21=_ϕ2-_ϕ1

sgnis the sign function, which delivers the sign of the value it receives.

It should be noted that (A1) and (A2) are dependent on each other, thus it requires a numerical method or software for its solution; in this case, the fsolve function of MATLAB ® was used.

Step 4. One of the reasons this method is widely used for mesh independence is because it calculates approximate values extrapolated from ϕ, to have a notion of the behavior of the variable as the number of cells increases or decreases. For this, the following equations are used:

_ϕext21=(_{r21p ϕ1}-_ϕ2)/(_r21p-1)

_ϕext32=(_{r32p ϕ2}-_ϕ3)/(_r32p-1)

Step 5. As a final step, the different errors are calculated and analyzed, as well as the mesh convergence index:

Relativeerror:_ea21=|_ϕ1-_ϕ2ϕ1|

Extrapolatederror:_eext21=|_ϕext21-_ϕ1ϕext21|

GC_I21=1.25_ea21r21p-1

This analysis is not performed for 32 because it is essential to analyze the refinement of the mesh, not the thickening of the mesh. In conclusion, according to NASA's Glenn Research Center, the GCI is a percentage measure of how far away the result of a simulation is from the asymptotic trend between meshes; thus, the smaller it is, the more constant the response will be and the better the discretization will be refined.

Processes 08 00878 g0a1 550

Figure A1.Isometric view of inclined: Rushton impeller (A); Paddles impeller (B); Propeller impeller: default (C); Propeller impeller: small design (D); and Brunswick Bioflo/CelliGen 115 bioreactor vessel (E).

Figure A1.Isometric view of inclined: Rushton impeller (A); Paddles impeller (B); Propeller impeller: default (C); Propeller impeller: small design (D); and Brunswick Bioflo/CelliGen 115 bioreactor vessel (E).

Processes 08 00878 g0a2 550

Figure A2.Glycerin viscosity as a function of the percentage of water.

Figure A2.Glycerin viscosity as a function of the percentage of water.

Processes 08 00878 g0a3 550

Figure A3.Statistical analysis: (top left) the main effects of dissolved oxygen, where the velocity label corresponds with agitation velocity [RPM] and the airflow has units of [L/min]; (top right) interaction effects for oxygen dissolved; (bottom left) data homoscedasticity; and (bottom right) Pareto chart of the standardized effects.

Figure A3.Statistical analysis: (top left) the main effects of dissolved oxygen, where the velocity label corresponds with agitation velocity [RPM] and the airflow has units of [L/min]; (top right) interaction effects for oxygen dissolved; (bottom left) data homoscedasticity; and (bottom right) Pareto chart of the standardized effects.

Processes 08 00878 g0a4 550

Figure A4._kLa determination for one agitation configuration: (A) 2.5lmi^n-1 and 100 RPM; (B) 2.5lmi^n-1 and 250 RPM; (C) 5lmi^n-1 and 100 RPM; and (D) 5lmi^n-1and 250 RPM.

Figure A4._kLa determination for one agitation configuration: (A) 2.5lmi^n-1 and 100 RPM; (B) 2.5lmi^n-1 and 250 RPM; (C) 5lmi^n-1 and 100 RPM; and (D) 5lmi^n-1and 250 RPM.

Processes 08 00878 g0a5 550

Figure A5._kLa determination for two impeller configurations: (A) 2.5lmi^n-1 and 100 RPM; (B) 2.5lmi^n-1 and 250 RPM; (C) 5lmi^n-1 and 100 RPM; and (D) 5lmi^n-1and 250 RPM.

Figure A5._kLa determination for two impeller configurations: (A) 2.5lmi^n-1 and 100 RPM; (B) 2.5lmi^n-1 and 250 RPM; (C) 5lmi^n-1 and 100 RPM; and (D) 5lmi^n-1and 250 RPM.

Table

Table A1.Experimental design.

Table A1.Experimental design.

StdOrder	RunOrder	PtType	Blocks	Impeller	Airflow	Velocity	DO (%)
12	1	1	1	Paddles	5.0	250	90.3
2	2	1	1	Propeller	2.5	250	52.2
3	3	1	1	Propeller	5.0	100	64.0
6	4	1	1	Small Propeller	2.5	250	52.2
13	5	1	1	Rushton	2.5	100	57.4
7	6	1	1	Small Propeller	5.0	100	43.5
9	7	1	1	Paddles	2.5	100	46.9
11	8	1	1	Paddles	5.0	100	64.9
5	9	1	1	Small Propeller	2.5	100	28.8
1	10	1	1	Propeller	2.5	100	47.6
16	11	1	1	Rushton	5.0	250	88.2
14	12	1	1	Rushton	2.5	250	84.9
15	13	1	1	Rushton	5.0	100	66.2
4	14	1	1	Propeller	5.0	250	70.3
10	15	1	1	Paddles	2.5	250	84.7
8	16	1	1	Small Propeller	5.0	250	66.8

Table

Table A2.Analysis of variance.

Table A2.Analysis of variance.

Source	DF	Adj SS	Adj MS	F-Value	p-Value
Model	12	4694.14	391.18	38.96	0.006
Linear	5	4234.85	846.97	84.35	0.002
Impeller	3	1803.46	601.15	59.87	0.004
Airflow	1	618.77	618.77	61.63	0.004
Velocity	1	1812.63	1812.63	180.53	0.001
2-Way Interactions	7	459.28	65.61	6.53	0.076
Impeller*Air flow	3	69.26	23.09	2.30	0.256
Impeller*velocity	3	373.42	124.47	12.40	0.034
Air flow*velocity	1	16.61	16.61	1.65	0.289
Error	3	30.12	10.04
Total	15	4724.26

Table

Table A3.Cells reference information.

Table A3.Cells reference information.

Parameter	Viscosity (Pa·s)	Max Shear Stress (Pa)	Max Shear Rate (s-1)	Cell Concentration(cellscm3)
EC	1.18×10-3	1292	1.09×106	2.81×108
SC	1.31×10-3	1250	9.54×105	1.20×109
CHO	2.80×10-1	0.10	3.58×10-1	4.44×105
TN-368	1.13×10-3	0.29	2.56×102	2×1011
SF-9	1.61×10-3	13.1	8.11×103	2×1011
AG	6.0×10-2	1.5	2.5×101	1×108
HeLa	1.0×10-3	2	2,0×103	6×105

References

1. Pérez, J.C.B.; Sánchez, R.A.H. Efecto de la relación agitación-aireación sobre el crecimiento celular y la producción de Azadiractina en cultivos celulares de Azadirachta indica A. Juss. Rev. Fac. Nac. Agron. 2010, 63, 5293-5305.

2. Post, T. Understand the Real World of Mixing. AIChE J. FG-1 2010, 106, 25-32.

3. Doran, P.M. Bioprocess Engineering Principles; Academic Press: Cambridge, MA, USA, 1995; ISBN 9780080528120.

4. Ochoa, F.; Gomez, E. García- & Bioreactor scale-up and oxygen transfer rate in microbial processes: An overview. Biotechnol. Adv. 2009, 27, 103-226.

5. Gogate, P.R.; Beenackers, A.A.C.M.; Pandit, A.B. Multiple-impeller systems with a special emphasis on bioreactors: A critical review. Biochem. Eng. J. 2000, 6, 109-144.

6. Buffo, M.M.; Corrêa, L.J.; Esperança, M.N.; Cruz, A.J.G.; Farinas, C.S.; Badino, A.C. Influence of dual-impeller type and configuration on oxygen transfer, power consumption, and shear rate in a stirred tank bioreactor. Biochem. Eng. J. 2016, 114, 130-139.

7. Uribe, A.R.; Rivera, R.; Aguilera, A.F.; Murrieta, E. Stirring and mixing. Enlace Químico 2012, 4, 22-30.

8. Yunus, A.; Cimbala, J.M.; Sknarina, S.F. Mecánica de Fluidos: Fundamentos y Aplicaciones; Primera, Ed.; McGrawHill: New York, NY, USA, 2001; pp. 10-11.

9. Centeno Carrillo, J.; Maroto Centeno, J. Utilización de un frasco de Mariotte para el estudio experimental de la transición de régimen laminar a turbulento. Rev. Española Física 1999, 13, 42-47.

10. Bird, R.B.; Stewart, W.E.; Lightfoot, E.N. Transport Phenomena, 2nd ed.; John Wiley & Sons: Hoboken, NJ, USA, 2006; ISBN 0470115394.

11. Deen, W.M. Transport in turbulence flow. In Analysis of Transport Phenomena; Oxford University Press: New York, NY, USA, 1998; ISBN 0195084942.

12. Chisti, Y.; Moo-Young, M. On the calculation of shear rate and apparent viscosity in airlift and bubble column bioreactors. Biotechnol. Bioeng. 1989, 34, 1391-1392.

13. Alfonsi, G. Reynolds-Averaged Navier-Stokes Equations for Turbulence Modeling. Appl. Mech. Rev. 2009, 62, 040802.

14. Proceedings of the STAR South East Asian Conference Best Practices Workshop: Heat Transfer, Singapore, 8-9 June 2015; p. 1.

15. Hanjalić, K.; Popovac, M.; Hadžiabdić, M. A robust near-wall elliptic-relaxation eddy-viscosity turbulence model for CFD. Int. J. Heat Fluid Flow 2004, 25, 1047-1051.

16. Domen, J.; Jain, R.; Cui, Z. Chapter 4-Environmental Impacts of Mining. In Environmental Impact of Mining and Mineral Processing; Butterworth Heinemann: Oxford, UK, 2016; pp. 1-53.

17. Collett, R.S.; Oduyemi, K. Air quality modelling: A technical review of mathematical approaches. Meteorol. Appl. 1997, 4, 235-246.

18. SIEMENS Using the Volume Of Fluid (VOF) Multiphase Model. Available online: http://mdx2.plm.automation.siemens.com/sites/default/files/Presentation/18% (accessed on 30 May 2020).

19. Armenante, P.M.; Nagamine, E. Uehara Effect of Low Off-Bottom Impeller Clearance on the Minimum Agitation Speed for Complete Suspension in Stirred Tanks. Chem. Eng. Sci. 1998, 9, 1757-1775.

20. Eppendorf. AG New Brunswick BioFlo ® /CelliGen ® 115 Operating Manual; Eppendorf: Hamburg, Germany, 2012.

21. Spanjers, H.; Olsson, G. Modelling of the dissolved oxygen probe response in the improvement of the performance of a continuous respiration meter. Water Res. 1992, 26, 945-954.

22. LabX New Brunswick Bioflo Reactors and Fermentors. Available online: https://www.labx.com/product/new-brunswick-bioflo-reactors-and-fermentors (accessed on 30 May 2020).

23. Celik, I.; Chen, C.J.; Roache, P.J.; Scheurer, G.; Washington, D.C. Quantification of Uncertainty in Computational Fluid Dynamics. Eng. Div. Summer Meet. 1993.

24. Eça, L.; Hoekstra, M. A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies. J. Comput. Phys. 2014, 262, 104-130.

25. Furukawa, H.; Kato, Y.; Inoue, Y.; Kato, T.; Tada, Y.; Hashimoto, S. Correlation of power consumption for several kinds of mixing impellers. Int. J. Chem. Eng. 2012, 2012, 1-6.

26. Spogis, N. Metodologia para determinação de curvas de potencia e fluxos caracteristicos para impelidores axiais, radiais e tangenciais utilizando a fluidodinamica computacional. Dissertation (master's), Universidade Estadual de Campinas, Faculdade de Engenharia Quimica, Campinas, SP. 2002. Available online: http://www.repositorio.unicamp.br/handle/REPOSIP/266231 (accessed on 5 June 2020).

27. Navisa, J.; Sravya, T.; Swetha, M.; Venkatesan, M. Effect of bubble size on aeration process. Asian J. Sci. Res. 2014, 7, 482-487.

28. Mcginnis, D.F.; Little, J.C. Predicting diffused-bubble oxygen transfer rate using the discrete-bubble model. Water Res. 2002, 36, 4627-4635.

29. Lange, H.; Taillandier, P.; Riba, J.-P. Effect of high shear stress on microbial viability. J. Chem. Technol. Biotechnol. 2001, 76, 501-505.

30. Keane, J.T.; Ryan, D.; Gray, P.P. Effect of shear stress on expression of a recombinant protein by Chinese hamster ovary cells. Biotechnol. Bioeng. 2003, 81, 211-220.

31. Iordan, A.; Duperray, A.; Verdier, C. Fractal approach to the rheology of concentrated cell suspensions. Phys. Rev. E 2008, 77, 11911.

32. Cai, M.; Zhang, Y.; Hu, W.; Shen, W.; Yu, Z.; Zhou, W.; Jiang, T.; Zhou, X.; Zhang, Y. Genetically shaping morphology of the filamentous fungus Aspergillus glaucus for production of antitumor polyketide aspergiolide A. Microb. Cell Fact. 2014, 13, 73.

33. Goldblum, S.; Bae, Y.K.; Hink, W.F.; Chalmers, J. Protective effect of methylcellulose and other polymers on insect cells subjected to laminar shear stress. Biotechnol. Prog. 1990, 6, 383-390.

34. Das, J.; Maji, S.; Agarwal, T.; Chakraborty, S.; Maiti, T.K. Hemodynamic shear stress induces protective autophagy in HeLa cells through lipid raft-mediated mechanotransduction. Clin. Exp. Metastasis 2018, 35, 135-148.

35. Chakraborty, S.; Nandi, S.; Bhattacharyya, K.; Mukherjee, S. Probing Viscosity of Co-Polymer Hydrogel and HeLa Cell Using Fluorescent Gold Nanoclusters: Fluorescence Correlation Spectroscopy and Anisotropy Decay. ChemPhysChem. 2020, 21, 406-414.

36. Karpinska Portela, A.M.; Bridgeman, J. Towards a robust CFD model for aeration tanks for sewage treatment-A lab-scale study. Eng. Appl. Comput. Fluid Mech. 2017, 11, 371-395.

37. Karpinska, A.M.; Bridgeman, J. CFD as a tool to optimize aeration tank design and operation. J. Environ. Eng. 2018, 144, 5017008.