Computational Fluid Dynamics for Advanced Characterisation of Bioreactors Used in the Biopharmaceutical Industry – Part II: Case Studies

Abstract

The first part of this series on characterisation of bioreactors in the biopharmaceutical industry using computational fluid dynamics presented a literature review to illustrate how characterisation can be performed and which process engineering parameters can be determined using computational fluid dynamics (CFD). In addition, experimental validation methods were presented, and an overview of typical hardware and software was also provided. In this second part, a selection of the authors’ research results will be used to demonstrate how the process characterisation of mechanically driven bioreactors for the biopharmaceutical industry can be determined with CFD and then experimentally validated. Three stirred tank bioreactors with different filling volumes and stirrers were used to demonstrate power input and oxygen transfer in single- and two-phase simulations. For wave-mixed and orbitally shaken systems, the fluid flow was transiently simulated and experimentally validated. In addition, the power input was also determined for both systems.

1. Introduction

Part one of this series showed that process engineering characterisation is essential for understanding and optimising bioreactors and bioprocesses for the biopharmaceutical industry. Process engineering characterisation usually includes determining the specific power input $P/V$, hydrodynamic stress $\tau$, mixing time ${\mathrm{\Theta }}_m$, volumetric oxygen mass transfer coefficient $k_{\mathrm{L}}a$ and, if applicable, the $N_{\mathrm{S}1}$ criterion. Determining these parameters by computational fluid dynamics (CFD) and the experimental investigation of the parameters, which serves as validation, were both described in part one. Furthermore, the literature review of $50$ publications showed that Ansys Fluent and OpenFOAM are among the most widely used software solutions for CFD applications in biotechnology. Using selected, practice-relevant examples from the authors’ research, the second part of the series now demonstrates how process engineering characterisation can be carried out using CFD and then be validated with experimental investigations. The two most frequently used software solutions were selected for the investigations. In Section 2, the specific power input is determined for two stirred reactors with different working volumes and stirrer configurations. This was then validated experimentally by measuring the electric current and the torque, which were then used to determine the specific power input. The determination of the $k_{\mathrm{L}}a$ value by coupling population balance modelling and CFD is also presented for a stirred bioreactor and compared with the gassing-out method. In Section 3, the Flexsafe RM 2 L basic single-use bioreactor bag from Sartorius AG is used to demonstrate how the specific power input and the hydrodynamic stress in a wave-mixed bioreactor system can be determined. Particle image velocimetry (PIV) was used for validation. In Section 4, the influence of the contact angle, which is material-dependent, on the specific power input in shake flasks is investigated. For this example, the sickle height was used for validation purposes and compared with experiments and semiempirical formulas.

2. Stirred bioreactor

In this case study, some of the methodological procedures described in part one of this series are used to numerically simulate three different stirred bioreactors to determine the specific power input and the $k_{\mathrm{L}}a$ value. The specific power input was determined for two different stirrers (stirrer diameter of 20 mm and 40 mm) according to Bauer et al. [1] for a bearing-free, magnetically driven 2 L benchtop system described by Schirmer et al. [2], based on Levitronix’s levitating impeller technology using the BPS-i30 drive (Figure 1) to validate the simulations described below. For this purpose, the torque was calculated using the known motor constant $K_{\mathrm{t}}$ with $1.13\mathrm{N}\mathrm{cm}{\mathrm{A}}^{-1}$ using eqs (1) and (2) and the required current for stirring:

Figure 1.

Computer-aided design (CAD) geometry of the setup and design of the bearing-free, magnetically driven 2 L benchtop system based on Levitronix’s levitating impeller technology (left) and numerically derived fluid flow field with a stirrer diameter of 40 mm (right).

$$M=K_{\mathrm{t}}\cdot I$$E1

$$P/V=\frac{2\cdot \pi \cdot n\cdot M}{V}$$E2

To determine the specific power input using a numerical simulation, the power input of the stirrer was determined based on the predicted fluid flow (see Figure 1) and the acting torque. For this purpose, Ansys Fluent was used with the realisable $k$-$\epsilon$-model [3]. Both the vessel walls and the impeller were treated as nonslip boundaries with standard wall functions and the axial velocity at the fluid surface was set to zero. The stirrer rotation was implemented using the multiple reference frame (MRF) method. Discretisation was performed using the first-order upwind scheme and pressure-velocity coupling using the COUPLED algorithm. The fluid domain was defined by an unstructured mesh of approximately $8\cdot {10}^6$ tetrahedrons. The specific power input for the working volumes of 2 L and 4 L was determined for the bacterial version of the Minifors 2 (Infors AG) with a total volume of 6 L, using OpenFOAM’s simpleFoam algorithm, the Gauss upwind scheme, the MRF method, and the $k$-$\epsilon$-model of Launder and Spalding [4]. Discretisation of the fluid domain was performed with $1.7\cdot {10}^6$ cells (2 L) and $3\cdot {10}^6$ cells (4 L), respectively. Validation of the numerical simulations was performed by determining the specific power input at the same working volumes using a torque meter to which the shaft and the stirrers were attached [1, 5, 6].

The results illustrated in Figure 2 show good correlation between the experimentally and numerically determined specific power inputs for both the magnetically driven system with the BPS-i30 drive and the Minifors 2. With the exception of the 40-mm diameter stirrer in combination with the BPS-i30, all the other configurations resulted in values $>5\mathrm{kW}{\mathrm{m}}^{-3}$, making them suitable for microbial applications.

Figure 2.

Double logarithmic representation of the numerically and experimentally determined specific power inputs of the magnetically driven system using the BPS-i30 drive with stirrer diameters of 20 mm and 40 mm as well as of the bacterial version of the Minifors 2 with working volumes of 2 L and 4 L.

In contrast to the previous ones shown, which were single phase and steady state, the oxygen transfer in another 2 L stirred bioreactor (modified HyPerforma glass bioreactor from Thermo Fisher Scientific Inc.) was determined transiently via a two-phase Euler-Euler simulation [7] (Figure 3 left). Besides the rotation of the stirrer which was realised via MRF, the aeration was defined by a gas inlet boundary condition at the sparger (a gas outlet is defined at the bioreactor lid). The liquid side mass transfer coefficient $k_{\mathrm{L}}$ value can be calculated according to eq. (3) as a function of the energy dissipation rate (Brüning’s [8] adapted formulation was used for these studies). The specific interface $a$ can be calculated according to eq. (4) [9]. Thus, the volumetric oxygen mass transfer coefficient $k_{\mathrm{L}}a$ can be determined, by the product of $k_{\mathrm{L}}$ and $a$. The class method was used to model the gas bubbles and their size distribution with the population balance modelling approach in this example. 24 gas bubble size classes were used, which according to our own investigations and literature values proved to be a good compromise between accuracy and computation time (the computation time increases exponentially as the class number increases) [7, 10, 11]. The model of Laakkonen et al. [12] with binary bubble breakup was used to model bubble breakup, and the model of Coulaloglou and Tavlarides [13] was used to model the coalescence behaviour. A review of the influence of different coalescence and breakup models can be found in Seidel and Eibl [7]. The interfacial models also have a significant influence on the predicted specific surface area. Different models were investigated by Seidel and Eibl [7], with the model of Schiller and Naumann [14] proving to be the best for calculating the drag force, the model of Tomiyama et al. [15] for calculating the lift force, and the model of Lamb [16] for calculating the virtual mass force. Using this model, a $k_{\mathrm{L}}a$ value of $\left(14.3\pm 0.6\right){\mathrm{h}}^{-1}$ was determined for a stirrer speed of 600 rpm and an aeration rate of 0.5 vvm (Figure 3 right). This result was verified experimentally using the gassing-out method according to the DECHEMA e.V. Working Group for Single-Use Technology [1, 17]. An eightfold determination resulted in a measured value of $\left(11.1\pm 0.2\right){\mathrm{h}}^{-1}$. This deviation can be explained not only by the discretisation error but also much more by the modelling error. The use of various models, some of which are based on empirical approximations, can lead to deviations from the experimental investigations. Only the MRF approach was chosen for the stirrer motion, which is less accurate than the dynamic mesh approach. Nevertheless, this method is suitable for estimating the $k_{\mathrm{L}}a$ value in stirred systems:

Figure 3.

$k_{\mathrm{L}}a$ value determination by means of CFD coupled with population balance modelling. CFD simulation of the aerated and stirred 2 L HyPerforma glass bioreactor (left). The colour map corresponds to the logarithmised relative gas fraction. The gas flow is visualised using line integral convolution. Temporal variability of the $k_{\mathrm{L}}a$ value determined by CFD simulation and its validation (right). Since the gassing-out method was used for the validation, no temporal statement about the $k_{\mathrm{L}}a$ value is available.

$$k_{\mathrm{L}}=\frac{2}{\sqrt{\pi }}\mathit{\epsilon \nu }{\left(\frac{D_{O_2}}{\nu }\right)}^{\frac{1}{4}}$$E3

$$a=\frac{A_{\mathrm{g}}}{V}=\frac{6\alpha }{d_{32}}$$E4

3. Wave-mixed bioreactor

The rocker-type motion of a wave-mixed bioreactor can be described by a harmonic oscillation function [18]. This allows the deflection angle $\phi$ at time point $t$ to be predicted by eq. (5), whereby ${\phi }_{\mathit{max}}$ corresponds to the maximum deflection angle and $\overline{\omega }$ to the angular velocity, which is calculated using the frequency $f$ (eq. 6). The frequency itself corresponds to the set rocking rate:

$${\phi }_t={\phi }_{\mathit{max}}\cdot sin\left(\overline{\omega }t\right)$$E5

$$\overline{\omega }=2\pi \cdot f$$E6

The motion can be realised with the help of a changing gravity vector $\overrightarrow{g}$ ($9.81\mathrm{m}{\mathrm{s}}^{-2}$) that is decomposed into its $x$ and $y$ components for a given deflection angle (Figure 4b). Thus, the gravity vector $g_x$ and $g_y$ can be calculated using eqs (7) and (8). Choosing a maximum deflection angle of ${10}^{\circ }$ and a rocking rate of 25 rpm leads to the outcomes in (Figure 4c and d):

Figure 4.

Representation of a) the 3D-scanned CAD geometry of the Flexsafe RM 2 L basic, b) the decomposition of the gravity vector into its $x$ and $y$ components at a deflection of $\phi$, c) and d) the decomposition of the gravity vector into its $x$ and $y$ components over time at a maximum deflection of 10${}^{\circ }$ and 25 rpm.

$${\overrightarrow{g}}_x=\overrightarrow{g}\cdot sin{\phi }_t$$E7

$${\overrightarrow{g}}_y=\overrightarrow{g}\cdot cos{\phi }_t$$E8

Another way to mimic the movement of the bag is to rotate or translate the mesh with a user-defined function according to eq. (9):

$${\overline{\omega }}_t=\overline{\omega }\left({\phi }_{\mathit{max}}\cdot \frac{\pi }{180}\right)\cdot cos\left(\overline{\omega }t\right)$$E9

Both of the above-mentioned methods for implementing the motion in the wave-mixed bioreactor were used for simulations of a two-phase air and water system in a Flexsafe RM 2 L basic single-use bioreactor bag from Sartorius AG (see used CAD geometry in Figure 4a). The simulation was performed transiently for liquid volumes of 0.5 L, 1 L, and 1.5 L in Ansys Fluent using velocity-pressure coupling, the SIMPLE algorithm, and the volume of fluid (VOF) model of Hirt and Nichols [19] with two Eulerian phases and a mesh consisting of $1.5\cdot {10}^6$ tetrahedrons. A no-slip assumption was chosen for the walls of the Flexsafe RM. Turbulence modelling was performed using the $k$-$\omega$ model. To achieve sufficient convergence for this procedure, time steps of ${10}^{-4}\mathrm{s}$ were chosen at the beginning of the simulation.

The CFD simulations were validated by comparing experimentally determined liquid level (Figures 5 and 6) and PIV measurements (Figure 7) to the simulations. This showed good qualitative and quantitative congruence of the geometric expression of the wave along the bag cross section for different volumes, deflection angles, and rocking rates. The relative deviations remained in the order of magnitude of less than 10% for all of the simulations, which can be attributed to a slight temporal offset of the images being compared. This resulted from the fact that, due to the large amounts of data generated during the numerical simulations, the simulation results were only saved every hundredth of a second. As a result, minimal deviations can occur in the comparisons of the wave characteristics at each full degree, because results might not always be available for the exact time point. Another cause of slight deviations is the use of a rigid, 3D-scanned bag shape for the simulations (Figure 4a), since in reality slight changes in shape can be observed and, depending on how the bag is fixed, small folds may also form, neither of which were taken into account in the simulation [20].

Figure 5.

Qualitative comparison of the fluid surfaces from the fluid level measurements (top) and the numerical simulations (bottom) at a) a fluid volume of 0.5 L, a maximum deflection angle of 6${}^{\circ }$ at 35 rpm, b) a fluid volume of 1 L, a maximum deflection angle of 6${}^{\circ }$ at 25 rpm, and c) a fluid volume of 1.5 L, a maximum deflection angle of 10${}^{\circ }$ at 25 rpm; all images shown have an instantaneous deflection of 3${}^{\circ }$.

Figure 6.

Dimensionless comparison of the fluid surfaces from the level measurements and the numerical simulations with different motion implementations for a fluid volume of 1 L, a maximum deflection angle of 6${}^{\circ }$ and an instantaneous deflection of 4${}^{\circ }$ at 25 rpm; simulation 1 corresponds to the motion implementation by means of decomposition of the gravity vector into its $x$ and $y$ components; simulation 2 was realised by means of the mesh motion.

Figure 7.

Qualitative comparison of the fluid velocities from the PIV measurements (top right) and the results of the CFD simulation (top left) for 1 L at 10${}^{\circ }$ and 25 rpm and a quantitative comparison (bottom) of the velocities along the introduced line with a distance of 10 mm above the ground at an instantaneous deflection of 5${}^{\circ }$; the respective legends indicate the fluid velocities from 0.0 m s${}^{-1}$ to 0.6 m s${}^{-1}$.

When comparing the results of the different motion implementations, it becomes clear that only marginal differences occurred between the two simulations (Figure 6). The relative deviations of the simulations to each other amounted to a maximum of 5%, which can also be attributed to the previously mentioned reasons as well as to the different mesh handling of Ansys Fluent, since, for example, the change of the mesh position in the Cartesian coordinate system must be taken into account [21]. Nevertheless, this way of handling the mesh appears to result in significantly shorter computation time when opting for the mesh motion option. Computation times of at least 12 weeks were required to reach quasi-stationary periods for gravity vector decomposition, while the computation time for mesh motion was only 4 weeks. Therefore, the latter should be considered the method of choice for future simulations. Validation by means of PIV measurements at specific deflection angles also proved successful. The planes from the PIV measurements were compared to the same planes from the simulations, and the velocity profiles along introduced lines were quantitatively compared. Consequently, good agreement can be assumed to indicate that the simulation results are of high quality. Based on the results, the periodically changing power input and the mean volume-averaged shear stresses could be determined (example shown in Figure 8), which with values of up to $250\mathrm{W}{\mathrm{m}}^{-3}$ and a maximum of $0.026\mathrm{N}{\mathrm{m}}^{-2}$ indicate good conditions for the cultivation of mammalian cells. For this reason, it comes as no surprise that they are most commonly used for mammalian cell inoculum production processes in feeding and, more recently, perfusion mode [22, 23, 24, 25].

Figure 8.

Characterisation of the specific power input of the numerical simulations at 1 L, 10${}^{\circ }$, and 25 rpm over an entire period (top) as well as representation of the maximum values of the mean volume-averaged shear stresses over an interval of one period (bottom).

4. Orbitally shaken bioreactor

The process engineering and numerical characterisation of orbitally shaken systems poses some challenges for the user. First, the system must be defined. There are different designs of the most common system, the shake flask, be it in shape (Erlenmeyer, Fernbach, Thomson, etc.), sensor installation, diameter of the flask, the neck, or in the material (various plastics and borosilicate glass) [26, 27, 28, 29]. In addition to the obvious parameters such as filling height and shaking rate, often neglected aspects such as the shaking amplitude or the contact angle of the material, which depends on the material properties (hydrophilic or hydrophobic), have a decisive influence on the oxygen transfer rate and the power input [30]. This can be explained by the fact that the power is transferred through the vessel wall into the suspension. Likewise, the liquid film that forms on the vessel wall plays a decisive role in the oxygen transfer, as this significantly increases the gas-liquid interface [31]. Finally, in orbitally shaken systems, effects such as the “out of phase” phenomenon can occur, whereas unfavourable cultivation parameters can cause an increase in shaking rate and lead to a decrease in power input, mixing efficiency, and oxygen transfer [32]. Therefore, great care must be taken when choosing simulation parameters and validating simulated orbitally shaken systems. In the first example, the numerically determined fluid motion and the height of the forming sickle are compared with empirical calculations and experiments (Figure 9). The simulations were performed using OpenFOAM v9, the VOF method with the interFOAM solver, and the k-$\omega$-SST turbulence model [34]. The CAD model of the $500\mathrm{mL}$ shake flask was created according to the geometry of Corning Inc.. For glass flasks, ISO 1773 and ISO 4797 could be used [35, 36]. The mesh created by snappyHexMesh consisted of 742,680 mesh cells. The simulation was run transiently for 20 s, using mesh motion (rather than manipulating the force vectors), with a maximum Courant-Friedrichs–Lewy (CFL) number of $0.9$. For experimental validation purposes, a flask filled with phenol red coloured water was shaken under identical conditions and the liquid movement was recorded. The calculation was based on an empirical estimation according to Büchs et al. [33], who developed a system of $16$ equations to be solved iteratively, that allow the maximal fluid height to be derived. The comparison of the maximum sickle heights showed very good agreement. This was 0.0461 m for the numerical approach, 0.0475 m for the experimental verification, and 0.0469 m for the empirical calculation.

Figure 9.

Simulation (left) and experiment (right) of the fluid flow in a 500 mL shake flask with 150 mL of working volume at 160 rpm and 25 mm shaking amplitude. The dashed red line indicates the empirically calculated maximum fluid height according to Büchs et al. [33], which deviates by only 2% from the simulated one.

Simulations of sickle height and general fluid motion are ideal for performing an initial verification of the accuracy of a simulation. The second example demonstrates the influence of the contact angle on the power input in shake flasks. The experimental determination of the power input in shake flasks has been described for a wide variety of parameter combinations by Büchs et al. [37, 38]. Glass flasks filled with water were primarily used, with contact angles between ${13}^{\circ }$ and ${44}^{\circ }$ being reported in the literature [39, 40]. Figure 10 compares experimentally measured data from 250 mL shake flasks (25 mL working volume and 25 mm shaker amplitude) to simulations with contact angles ranging from ${30}^{\circ }$ to ${75}^{\circ }$. Büchs et al. [37] determined specific power inputs from $0.1\mathrm{kW}{\mathrm{m}}^{-3}$ to $3\mathrm{kW}{\mathrm{m}}^{-3}$, while the results from the CFD simulations ranged from $0.04\mathrm{kW}{\mathrm{m}}^{-3}$ (${30}^{\circ }$, 120 rpm) to $4.4\mathrm{kW}{\mathrm{m}}^{-3}$ (${75}^{\circ }$, 360 rpm). Good agreement between simulation and measurement was shown at a contact angle of ${30}^{\circ }$. With higher contact angles, the power input increased, which must be taken into account when changing from glass to (single-use) plastic flasks.

Figure 10.

The specific power input for 250 mL shake flasks, which were shaken with an amplitude of 25 mm and different shaking rates. Different contact angles were investigated by CFD and compared with the experimental data from Büchs et al. [37].

5. Conclusions

The first part of the series illustrated how process engineering parameters for bioreactors used in the biopharmaceutical industry can be determined using CFD and how they can be validated. Furthermore, it was shown that CFD has become well established as a tool for characterising bioreactors, and that both hardware and software are constantly being improved to speed up the calculations and achieve higher accuracy. In this part of the series, practical application examples were used to demonstrate that with the appropriate choice of model, process engineering parameters can be determined to the required degree of accuracy (Table 1). This is independent of the mechanical drive for stirred, orbitally shaken, and wave-mixed bioreactor systems. Using different stirred bioreactors, it was possible to show that the calculated specific power input agrees with both the electrical and the torque measurement method. In addition, it was also demonstrated that the $k_{\mathrm{L}}a$ value in a stirred system can be calculated in a way that takes coalescence and bubble breakup into consideration, and can be validated with the gassing-out method. Both the free surface and flow velocities in a wave-mixed system were determined using CFD and validated experimentally. In addition, the specific power input and the hydrodynamic stress were also determined for this system. Using the example of shake flasks, it was shown that the free surface can also be precisely determined for orbitally shaken systems, with the sickle height of the CFD and experimental validation differing by less than 2%. Furthermore, the influence of the contact angle on the specific power input was also examined. As shown in this series, CFD is used extensively for the characterisation of bioreactors and is likely to be used even more in the future.

Bioreactor	Volume	Solver	Turbulence	Mesh in; ${10}^6$ cells	Boundary
OpenFOAM
Minifors 2	2 L${}^{\mathrm{a})}$; 4 L${}^{\mathrm{a})}$	simpleFoam${}^{\mathrm{d})}$	$k$-$\epsilon$	1.7${}^{\mathrm{f})}$; 3.0${}^{\mathrm{f})}$	Wall: no-slip; Surface: symmetry
HyPerforma	2 L${}^{\mathrm{b})}$	reactingTwo; PaseEulerFoam${}^{\mathrm{e})}$	$k$-$\epsilon$	0.3${}^{\mathrm{f})}$	Wall: no-slip and; Gas in−/outlet
Shake flask	0.025 L${}^{\mathrm{c})}$; 0.15 L${}^{\mathrm{c})}$	interFoam${}^{\mathrm{e})}$	$k$-$\omega$-$\mathit{SST}$	0.09${}^{\mathrm{f})}$; 0.7${}^{\mathrm{f})}$	Wall: contact angle
Fluent
BPS-i30	2 L${}^{\mathrm{a})}$	COUPLED${}^{\mathrm{d})}$	$k$-$\epsilon$	8.0${}^{\mathrm{g})}$	Wall: No-slip; Surface: symmetry
Flexsafe RM	0.5 L to; 1.5 L${}^{\mathrm{c})}$	SIMPLE${}^{\mathrm{e})}$	$k$-$\omega$	1.5${}^{\mathrm{g})}$	Wall: no-slip

Table 1.

Overview of the investigated bioreactor systems with working volumes, solver used, turbulence models, number, and type of mesh cells as well as boundary conditions. The systems are divided into OpenFOAM and fluent cases according to the simulation software used.

a) single-phase b) two-phase Euler-Euler c) VOF d) steady-state e) transient f) mostly hexahedra g) tetrahedrons

Acknowledgments

We would like to thank Sören Werner for contact angle investigations of the shake flasks and Darren Mace for English proofreading.

Conflict of interest

The authors declare no conflict of interest.