Oxygen Mass Transfer in Biopharmaceutical Processes: Numerical and Experimental Approaches

Oxygen supply in aerobic bioprocesses is of crucial importance. For this reason, this paper presents the oxygen demand of different cells and summarizes experimental and numerical possibilities for the determination of oxygen transfer in bioreactors. The focus lies on the volumetric oxygen mass transfer coefficient (kLa) calculation using computational fluid dynamics and state-of-the-art models for surface-aerated and forced-aerated bioreactors. In addition, experimental methods for the determination of the kLa value and the gas bubble size distribution are presented.


Introduction
Biopharmaceuticals are a fast-growing market with enormous potential for the pharmaceutical industry. The global market value is currently above $300 billion [1] and is forecasted to increase to $389 billion by 2024 [2]. Depending on the type of biopharmaceutical, either cell cultures of plant or animal origin or microorganisms are used. Since the majority of biopharmaceuticals are produced in aerobic processes, the oxygen supply of the cells plays a crucial role. Oxygen must continuously be present in dissolved form in the medium for a continuous utilization by the cells. Due to its low solubility, oxygen is the only nutrient that needs to be added throughout the whole process. If oxygen is limited or the oxygen concentrations are fluctuating, the growth rate is reduced and may lead to a reduced product formation rate and changes in metabolism [3][4][5]. The amount of oxygen required by the cells can be expressed by the specific oxygen uptake rate q O 2 , which is characteristic for each organism, and may vary in different process phases and under different process conditions [6].
The oxygen required is supplied to the respective bioreactor system in the form of air, pure oxygen or a mixture of both, whereby the type of aeration varies according to the system. Stirred tank bioreactors (STR), which are the most frequently used bioreactors in biotechnology [7], are typically gassed by forced aeration. In contrast, orbitally shaken systems and wave-mixed bioreactors are usually only gassed via the surface of the liquid.
To meet the oxygen demand of the cells, the oxygen transfer rate (OTR) must be equal to or greater than the oxygen uptake rate (OUR) of all cells, whereby the OTR is calculated from the product of k L a and the oxygen concentration difference. A variety of empirical formulas and correlations for calculating the OTR and k L a are described in the literature [8][9][10][11]. For the experimental determination of the k L a value, different methods have been developed [12,13]. The values determined can vary considerably, depending on the measurement method.
Furthermore, the surface area available for gas exchange can also be investigated. Therefore, the specific interfacial area is of particular interest for polydisperse systems such as STR and airlift bioreactors, since the gas exchange surface is proportional to the oxygen transfer rate.
To characterize bioreactors in terms of process engineering, computational fluid dynamics (CFD) simulations are increasingly used [14][15][16]. CFD allows time-and spaceresolved insights into bioreactors and to optimize them in silico already in the development phase. Using two-phase modeling, it is possible to make predictions about the volumetric oxygen mass transfer coefficient for all types of cultivation systems. The choice of different CFD models and the eventual coupling with population balance models (PBM) is of crucial importance for an economical and realistic prediction [17,18].
This review evaluates the potential of experimental and numerical approaches to the prediction of the oxygen transfer rate, with special emphasis on CFD-based methods. Firstly, an overview of specific oxygen uptake rates of different organisms used in the biopharmaceutical production is presented. This is followed by a summary of empirical calculation methods for the k L a value in different bioreactor systems. Finally, the possibility to numerically determine the k L a value by CFD will be described and discussed. Different models for the respective bioreactor systems will be presented, and the economic efficiency of the model selection and the simulation infrastructure will be evaluated. In addition, the influence of mass transfer, interface, coalescence and bubble breakup models on the calculation of k L a will be presented.

Oxygen Demand of Cells
To ensure a sufficient supply of oxygen to the cells (this term includes animal and plant cells as well as microbials), the OTR in the bioreactor must be greater than or equal to the OUR (Eq. (1)).

OTR OUR
(1) The OTR will be discussed in detail in Sect. 3. The OUR is the product of the cell-specific oxygen uptake rate q O 2 and the biomass concentration c x (Eq. (2)).
The cell-specific oxygen uptake rate is assumed to be constant during the exponential growth [19] phase but changes within other process phases such as lag or plateau phase [6]. Furthermore, the growth rate [20] and the selection of the limiting carbon sources [21] affect the specific oxygen demand. Thus, a reasonable value for q O 2 , depending on the cell line, the process and the media must first be determined or taken from suitable literature sources.

Determination of Specific Oxygen Uptake Rate
For the experimental determination of the specific oxygen uptake rate, several approaches exist, which have been well explained by [8,22,23]. The most-used techniques are summarized in Tab. 1.
Applying a gas phase mass balance calculation requires precise oxygen analyzers for the gas inlet and outlet as well as controlled flow rates, e.g., by mass flow controllers. This approach is suitable for all types of bioreactors but may be less accurate if the gas flow rate is high and the q O 2 and bioreactor volume are relatively low. The gas concentrations measured can be used to calculate the OTR and thus the OUR. To finally calculate the q O 2 , the OUR must be divided by the biomass concentration, determined for example by an offline measurement. The dynamic evaluation approach is based on the measurement of the oxygen consumption within the liquid phase. During the running process, the oxygen supply is switched off and the decrease of the dissolved oxygen concentration is recorded. It must be ensured that this measuring phase does not take too long in order to avoid oxygen limitation in the process. The slope of the dissolved oxygen decline can then be divided by the biomass concentration to calculate q O 2 . The measurement can be performed several times in the process. For the dissolved oxygen concentration profile approach, either the OTR must be known, or the k L a value and dissolved oxygen concentration in order to calculate it. With this, the OUR can be then calculated (Eq. (5)) and thus the q O 2 . However, it must be noted that the measured k L a values also differ depending on the measurement method, so that errors and inaccuracies in the k L a value determination affect the calculation of the q O 2 accordingly.
Further, different methods like the yield coefficient approach exist, which uses a simplified stochiometric calculation on the basis that the OUR equals the product of biomass concentration and growth rate divided by the cell yield on oxygen [8]. The accuracy of this method can be increased by introducing the maintenance coefficient, as the overall yield is not constant during the process. An overview of yield coefficients of different microorganisms is provided by Heijnen and Roels [280].
values determined by dissolved oxygen concentration and OTR [8].

Literature Values for the Specific Oxygen Uptake Rate
The specific oxygen uptake rate has already been investigated and analyzed by countless authors. Therefore, only selected summaries are presented in the following and the most important findings are discussed. A comprehensive overview of the oxygen demand of different cell lines is provided by [6]. Variations of dissolved oxygen concentration, temperature and pH for Chinese hamster ovary (CHO) cell cultivations were analyzed by [26]. An overview of oxygen uptake rates of Spodoptera frugiperda cells of subclone 9 (Sf9) in different media is provided by [67]. From the available literature sources, it can be seen that the specific oxygen uptake rates vary greatly due to the use of different measuring methods, clones, media and processes and, therefore, cannot always be adopted one-to-one for the process of interest. The distribution of the q O 2 values for different cell lines from the literature references considered is summarized in Tab. 2. In general, it is recommended to measure the oxygen demand of a bioprocess itself. As an alternative, the third quartile of the q O 2 value distribution can be taken from Tab. 2 and a corresponding k L a value calculated, with which the OTR would be sufficiently high in 75 % of all experiments considered.

Oxygen Mass Transfer
As oxygen is barely soluble in water and permanently consumed by the cells, oxygen must be added to the system continuously. Oxygen passes from the gas phase through several steps and resistances into the cells where the oxygen is consumed by cellular respiration [12]. The individual steps are described in detail in [8,82,83], typically using the two-film theory of Whitman [84], which is discussed in Sect. 3.1 in more detail. An overview of alternative gas transfer theories can be found in [85].

Two-Film Theory
The total coefficient of liquid phase mass transfer K L , from the two-film theory, which according to Eq. (3) is composed of the liquid side mass transfer coefficient and the gas side mass transfer coefficient, can be simplified. Because of the low solubility of oxygen in water k G is much larger than k L ; thus, the k G term can be neglected and K L k L [83,86,87]. The temperature dependence is described by linking the temperature with the universal gas constant and the Henry constant H cp in Eq. (3). If oxygen is considered to be an ideal gas, the dimensionless Henry constant H cc can also be used due to the connection of H cc = R T H cp . The OTR is often used to quantify the oxygen transfer. The OTR is the product of the liquid side mass transfer coefficient k L , interfacial area a and the oxygen concentration gradient (Eq. (4)). In experiments it is only possible to measure the total mass transfer K L and not separately the individual mass transfer coefficients k L and k G , but by the simplification described previously K L can be assumed as k L . Nevertheless, it is often challenging to measure k L and a individually in practice, therefore, the two parameters are combined to form k L a [8]. Measurement methods for k L a and a are described in Sect. 4.1 and 4.2. The dependence of OTR on hydrodynamic parameters is illustrated in [12]. The OTR depends on physical properties of fluids, bioreactor geometry, process parameters and biomass concentration. The change of oxygen concentration in the system can be described by Eq. (5). To avoid oxygen limitation during cultivation, the OTR must be greater than or equal to the OUR.

Bubble Breakup and Coalescence
In the case of forced gassing, the specific interfacial area a depends mainly on the size and shape of the gas bubbles. Thus, two processes have a decisive influence on the gas bubble size distribution and, hence, on the available mass transfer area. Firstly, gas bubbles can break up into two or more gas bubbles, and secondly, gas bubbles can also merge, processes called breakup and coalescence, respectively [95]. Four mechanisms that can lead to gas bubble breakup are described in the literature. In turbulent flows, such as in bioreactors, breakup due to turbulent fluctuations and collisions is most frequent. The other three mechanisms are viscous shear stress, shearing-off processes and interfacial instability [96]. The modeling of bubble breakup typically includes two submodels, one for the breakup frequency and one for the resulting daughter bubble size distribution. Most models assume a binary breakup [95].
The second phenomenon besides the bubble breakup is the more complicated coalescence [97]. In addition to the bubble-liquid interaction, bubble-bubble interaction is also involved in this process. In the literature, three different theories that describe the coalescence process are mentioned. The earliest and most frequently used theory is the film-drainage-model of Shinnar and Church [98]. According to this theory, there are three steps until complete coalescence occurs. First, two gas bubbles meet and trap a small amount of liquid in the form of a film. As long as the gas bubbles remain in contact, liquid drains from the film. If the liquid level falls below a critical thickness, the film breaks and the gas bubbles fuse. Coalescence only occurs if the interaction time is long enough to displace the critical amount of water from the film.
The second theory comes from Howarth [99]. This is the energy model, which assumes that if the collision energy is high enough, the gas bubbles coalesce directly without first being separated by a liquid film. The third and newest theory is from Lehr et al. [100]. This critical velocity model is an empirical theory based on experimental observations. For coalescence to occur at all, gas bubbles must first collide. Typically, a collision occurs if the relative velocities of the gas bubbles differ. According to Liao and Lucas [101], the occurrence of relative velocity differences can be divided into five different categories: turbulence fluctuation, viscous shear stress, capture in turbulent eddies, buoyancy and wake interaction. For the modeling of the coalescence frequency, either a fully empirical model is used, or one based on physics. In physical models, the coalescence frequency is derived from the product of coalescence efficiency and collision frequency.

Experimental Approaches
Different experimental methods are used to determine the oxygen transfer and the gas bubble size distribution. An overview of possible methods is described in this section. They can be used for the process engineering characterization of bioreactor systems or for the validation of CFD simulations (Sect. 5).

Volumetric Oxygen Mass Transfer Coefficient Measurement
Over the years, various methods for measuring the volumetric oxygen mass transfer coefficient (k L a) have been developed and extensively investigated, whereby a distinction can be made between abiotic and biotic measurement methods. Usually, bioreactors are first characterized abiotically (i.e., in the absence of cells) before a biological characterization is performed during cultivation. The following methods are the most common: -abiotic saturation-based methods (gassing-out method and dynamic sulfite method) [103,104], -measurement with chemical model system (static sulfite method, hydrazine method, glucose oxidase method) [13], -biotic saturation-based methods (respiratory gassing-out method) [102], -biotic steady-state methods (gas balancing method) [12].
The most common measurement method utilize nitrogen to strip out oxygen. After the oxygen is stripped out, an aeration with oxygen follows again. The slope of the saturation process, which occurs by the aeration with oxygen, can be used to determine the OTR. Since there are no cells active during saturation, OUR is zero and can be ignored. The method can be used with water or cell culture media in order to reach physical fluid properties as closely as possible. Similarly, the dynamic sulfite method utilizes sodium sulfite to deplete the oxygen and the saturation slope is measured. The reaction is catalyzed by a heavy metal ion, mostly copper, which must be in a very tight range. Furthermore, heavy metal ions need to be disposed of correctly. The methods described, based on oxygen saturation after an oxygen stripping by either cells, nitrogen or chemical depletion, are heavily dependent on the exact measurement of the dissolved oxygen concentration. Classical amperometric oxygen probes and also modern probes based on a membrane with a fluorescent dye exhibit a certain response time, which influences the measurement and even makes it impossible to measure in the range of high k L a values, since oxygen saturation is faster than the probe response time. Another limitation is the influence of the gas bubbles on the probe measurement during saturation, which is eliminated by a measurement within an external loop, eliminating all bubbles with a membrane. In addition, the specific location of the measurement influences the result and location change of measurement is complicated or impossible due to the bioreactor geometry and probe ports.
The respiratory gassing-out method is applied during an in vitro cultivation. The oxygen supply is switched off during the cell growth phase, and thus, oxygen is depleted. The dissolved oxygen concentration decreases at a rate equal to the OUR. Thus, the following simplified Eq. (6) can be used.
When switching on oxygen supply before reaching the critical oxygen concentration of the cell type used, the oxygen concentration increases until saturation is reached. The slope of the curve during the oxygen depletion can be used to determine the OUR. The slope of the following saturation process together with the determined OUR can then be used to calculate the OTR and thus, the k L a value [12,103,105]. In comparison to the respiratory gassing-out method, oxygen concentration in the gas balancing method is additionally measured for the inlet and outlet gas stream by a gas analyzer, resulting in an oxygen mass balance under steady-state conditions. At steady state, the OTR corresponds to the OUR. Therefore, k L a can be calculated according to Eq. (7).
However, the gas balancing method demands very sensitive measurement, in particular with smaller bioreactors, since the difference between incoming and output gases can be very small. Still, it is considered to be the most exact method, since the cells consume the oxygen while oxygen is supplied, and thus, the original state of the cultivation is used to determine the gas mass transfer.
It should be noted that the choice of measurement method has a considerable influence on the measured k L a values [12]. For example, measurements with the sulfite method achieved higher k L a values than those with electrolyte solution, especially at low stirrer speeds [12]. A direct comparison of the gassing-out method with the gas balancing method also showed significant differences, although no clear trend could be identified [102].

Bubble Size Distribution Measurement
In contrast to surface-aerated systems, the determination of the specific interphase area is much more complicated for polydisperse systems. The interphase area is determined experimentally through the gas bubble size distribution, whereby the measuring methods can be differentiated into intrusive and extrusive measuring methods. Extrusive methods offer the advantage that the measurement does not influence the flow, and often the whole bioreactor can be investigated, but these measurement methods often reach their limits at high aeration rates. Furthermore, the bioreactors examined must be optically transparent, which leads to problems, especially with stainless-steel systems. With intrusive methods, the measurement can often only be carried out locally. Various methods have been developed for solid particles, but certain methods can also be used for gas bubbles, although these are physically different from solid particles. A detailed overview of different methods for measuring particle size distributions can be found in [106,107].
Examples of intrusive methods include focused beam reflectance measurement (FBRM) [108][109][110] and photooptical methods such as capillary suction probe technique (CSP) [111][112][113] or the smart online particle analysis technology (SOPAT) [114,115]. The FBRM method is suitable for high bubble densities and measures the chord length distribution via the backscatter of a rotating infrared laser [116]. With the SOPAT probe, which is based on the endoscope principle, the two-dimensional image of the gas bubbles in the measuring window is recorded continuously [114]. In CSP, gas-liquid dispersion is extracted via a transparent capillary under vacuum. In the capillary, the gas bubbles are measured by a LED phototransistor pair [111].
Essential extrusive methods include laser doppler anemometry (LDA) [117][118][119], ultrasonic [120][121][122][123][124] and shadow imaging (shadowgraphy) [125][126][127][128]. For LDA measurement, two laser beams cross each other in the transparent bioreactor. If bubbles move through the beams, the light is scattered, and the particle velocity can be determined via the scattering. Per velocity coordinate one laser and one detector are required. Accordingly, three double beam lasers could be used to measure the three-dimensional bubble velocity in the STR. However, the bubble size is only indirectly determined by the correlation of bubble size and bubble rising velocity. This correlation is only reliable for vertically rising gas bubbles and, therefore, the gas bubble size determination in STR using LDA is not suitable. More fundamental details about LDA can be found in [129]. Ultrasonic-based methods typically consist of an ultrasonic transmitter and an ultrasonic receiver. The attenuation of the sound is used to measure the size of the gas bubbles. The ultrasonic method offers the great advantage that the bioreactor does not have to be transparent. It has been shown that comparable results can be achieved with ultrasonic and optical methods [121]. Shadowgraphy measurements belong to the optical measurement methods. Accordingly, the gas bubble size distribution can only be measured in transparent bioreactors. Like the LDA and ultrasonic method, the measurement setup consists of an emitting device and a receiver on the other side of the bioreactor. In the shadowgraphy method, diffuse light is emitted onto the bioreactor and the shadows created by the gas bubbles are recorded by a camera (Fig. 1). Typically, cylindrical bioreac-tors are placed in a rectangular transparent container and surrounded with water to minimize optical distortion [130]. Due to a high gassing rate, the number of bubbles exceeds the capacity for optical methods such as shadowgraphy, as bubble overlapping and swarms of bubbles occur [131]. With the help of sophisticated image processing algorithms, it is nevertheless possible to measure gas bubble size distributions in such systems [131]. For the image evaluation, classical image processing algorithms like filter and watershed algorithm [132,133] are used as well as more and more convolutional neural networks [134,135].

Numerical Approaches
Numerical fluid dynamics have been used for years for the process engineering characterization of bioreactors. Compared to physical experiments, CFD offers more detailed results, is relatively fast, reproducible and usually cheaper [18]. Furthermore, parameter ranges can be investigated with relatively little effort and are not dependent on prototypes. Numerical methods can be used to calculate the flow field, shear gradients, mixing times, power input and sedimentation criteria. CFD is based on the conservation laws of mass, momentum and energy. Fluid dynamic problems can be solved using the Navier-Stokes equations. Fundamental information on the theoretical background of CFD can be found in various textbooks [18,136,137].
In industrial applications, including biotechnology, turbulent flows are often required to achieve sufficient mixing.
Turbulent flows place special demands on numerical simulation and can be modeled with varying degrees of accuracy [18]. In direct numerical simulation, the conservation equations are solved numerically without turbulence modeling. In order to reduce the calculation effort, various turbulence models have been developed. In large eddy simulation, the conservation equations are solved only for large eddy structures directly. Models approximate smaller eddies. The Reynolds-averaged Navier-Stokes approach (RANS) is the most frequently used approach for turbulence modeling. The RANS approach does not solve turbulence structures at all, which makes it not only the simplest, but also the fastest calculation method. Therefore, in the following only RANS approaches will be considered. With the help of Figure 1. Schematic representation of the three extrusive methods shadowgraphy, LDA and ultrasonic. For the shadowgraphy method, the bioreactor is illuminated with diffuse light from one side, and the shadows of the gas bubbles are recorded on the opposite side by the chargecoupled device (CCD) camera and evaluated by image processing software. A laser in combination with a diffuser can be used as a light source. In the LDA method, the emitted laser beam is split using a beam splitter. The beams cross inside the bioreactor through a lens. The volume measured is located in this area of the bioreactor. The Doppler-shifted frequency can then be recorded by the photomultiplier tube and subsequently evaluated. For the ultrasonic method, no transparent bioreactor is necessary. By using a ring of transmitter and receiver, the horizontal bioreactor plane can be monitored, and multiple gas bubbles can be measured simultaneously.
The bubble size can be determined by evaluating the time of flight.
two-phase models, it is possible to determine further process engineering parameters, such as the volumetric oxygen mass transfer coefficient.

Two-Fluid Model
Over recent years, different modeling approaches have been developed for two-phase or multiphase modeling. Depending on the bioreactor type, certain models are more suitable for determining the k L a value than others. The two-phase approaches include the segment method, the marker-andcell method, level-set method, Monte-Carlo method and Lattice-Boltzmann method. In practice, the volume of fluid (VOF) method for continuous phases has become established. For disperse systems, the Euler-Lagrange and Euler-Euler approach are mainly used, whereby the latter is state of the art [107,138]. An overview of the state of the art of multiphase modeling can be found in [139].

Volume of Fluid
The VOF method was first described by Hirt and Nichols for free air-liquid interfaces and uses a mixed fluid [140,141]. The properties of the mixed fluid are weighted according to the phase fraction a i of the individual fluids (Eqs. (8) and (9)). Almost the whole system contains pure fluids; only at the phase interface does the mixed fluid have to be considered. The phase fraction transport equation is calculated using Eq. (10), whereby the condition Eq. (11) must be met. The continuity equation then corresponds to Eq. (12) and the momentum Eq. (13). F is the surface tension acting on the gas-liquid interface. For modeling the surface tension, several models exist, where besides the sharp surface force model and the parabolic reconstruction of surface tension model, especially the continuum surface force (CFS) and the smoothed CFS models are used [142][143][144][145]. The turbulence modeling is not considered here and can be found in [18,146]. Often variants of the k-e and k-w model are used, which both belong to the RANS models [147][148][149]. In both models two coupled transport equations are used to describe the turbulence. In the k-e model the turbulent kinetic energy k and the turbulent energy dissipation rate e are used. For the k-w model the specific rate of dissipation w is used instead of the turbulent energy dissipation rate. The VOF method is the most common method to model all kinds of surface aerated bioreactors. Shake flask, baffled and unbaffled 500 mL / 50-150 mL Ansys CFX 11.0 [147] Orbitally shaken bioreactor b) 600 mL / 100-500 mL Ansys Fluent 15.0 [148] Orbitally shaken bag 2 L / 1 L OpenFOAM 2.0.X [39] Conical bottom shaken vessel 6.3 L / 1-2.5 L Ansys Fluent 14.5 [155] Wave-mixed bioreactor e) 10 L / 5 L Ansys Fluent 13.0 [156] Orbitally shaken bag 20 L / 10 L OpenFOAM 2.0.X [39] Wave-mixed bioreactor e) 20 L / 2-10 L OpenFOAM (n.a.) [157] a) bioREACTOR (2mag); b) TubeSpin bioreactor (TPP); c) glass Erlenmeyer flask; d) plastic; e) cellbag (Cytiva, formerly GE Healthcare Life Sciences); n.a. not available suitable because fewer equations have to be solved compared to the Euler-Euler method and is therefore faster. In general, comparable k L a values to experimental data can be simulated [148][149][150]. Zhu et al. [148] showed an average deviation of 10.54 % relative to the static gassing-out method for different filling volumes and shaking rates.

Euler-Euler Approach
In contrast to the VOF method, the Euler-Euler method, which can be used for dispersed multiphase systems, does not use a mixed fluid, but the balance equations are solved individually for each fluid. The continuity equations for two fluids correspond to Eq. (14). However, there is no phase transition in this equation [158]. The momentum equation for two fluids is also shown without phase transition in Eq. (15) [158]. F I,i is the sum of all interfacial forces acting on phase i from the other phase [159]. These forces are not considered in the VOF method and include drag force, lift force, virtual mass force, wall lubrication force and turbulent dispersion force. Since it is often analytically impracticable to determine these forces exactly, partial models are used.

Drag Force
In the literature, the drag force is described as the interfacial force that has the most significant influence on the simulation [160,161]. The drag force is defined according to Eq. (16) [162,163]. The drag force coefficient C D is difficult to calculate analytically. For solid spherical bodies in laminar flow, this can be done according to Stokes' law [164]. Typically, therefore, models are used for the calculation of C D in practice. The most commonly used drag force coeffi-cient models are those proposed by Schiller & Naumann [165] and Ishii & Zuber [166,167]. Both models do not consider the influence of turbulence. Although this makes them attractive in terms of calculation time, it could also explain why these models often do not correspond to experiments [168,169]. The model of Gidaspow et al. [170] is a combination of the models of Ergun [171] and Wen & Yu [172] and covers a broader range of possible applications, whereby a higher computing time has to be accepted [173]. In contrast to the models by Ishii & Zuber [166] and Schiller & Naumann [165], this model takes the influence of turbulence into account. The model of Laín et al. [174] also takes turbulence into account, which makes the model attractive for STR. The model of Tomiyama et al. [175] is also more sophisticated than the standard models and considers not only the turbulence but also the aspect ratio. Gradov et al. [160], Basavarajappa & Miskovic [176], as well as Karimi et al. [169] were able to demonstrate that the model of Lane et al. [177] showed the best agreement with the experiments in turbulent flow. In addition, some models adjust the drag coefficient by swarm correction when particles are found in swarms [178]. A study on the influence of the drag model on the gas-liquid flow can be found in [179].

Lift Force
The lift force is another force that influences the momentum transfer. The lift force is defined by Eq. (17) and is the force acting on the dispersed particle perpendicular to the relative velocity direction [180]. Zhang et al. [181], who have carried out investigations on interfacial forces in stirred and aerated bioreactors, were able to show that the lift forces in the bioreactor are about 1000 times smaller than the drag force. This coincides with several other authors who also consider the lift force to be negligible [161]. Lucas et al. [182] were able to show that Tomiyama's model [183] was in good agreement with bubble column experiments. However, the lift force is still not fully understood and requires further numerical and experimental investigations [182].

Virtual Mass Force
The virtual mass force is a force that acts on the dispersed particle and has the same effect as an additional mass. The virtual mass force is defined according to Eq. (18). For spherical particles, the virtual mass force coefficient C VM is 0.5; therefore, many authors have assumed that C VM is 0.5 [163,184]. However, different models exist for this coefficient. The model of Pudasaini [185] is an up-and-coming model because it is a complete analytical model. The phenomenon of virtual mass is almost ubiquitous and, as Zhang et al. [181] have shown, is of great importance,  [186].

Wall Lubrication Force
The wall lubrication force was first taken into account by Antal et al. [187] (Eqs. (19) and (20)). It is the repulsive force of the wall, acting on a dispersed particle [188]. This effect is caused by an asymmetrical flow around the particle near the wall. Antal et al. [187] use the two constants C w1 = 0.01 and C w2 = 0.05 to calculate the wall lubrication coefficient. In addition to the model of Antal et al. [187], the models of Tomiyama et al. [189] and Frank [190] are widely used.

Turbulent Dispersion Force
The turbulent dispersion force was first considered by Lopez de Bertodano et al. [191] (Eq. (21)) and described the force that is transferred to the disperse phase by turbulent fluctuations in the continuous phase [182]. Lopez de Bertodano et al. [191] used values from 0.1 to 0.5 for the turbulent dispersion coefficient. Gosman et al. [192] established a relationship between the turbulent dispersion coefficient and the Stokes number. The latest model developed by Burns et al. [193] is based on a time-averaged resistance model. The turbulent dispersion force is often neglected in the literature. However, if the force is taken into account, the Gosman et al. model [192] is usually used because it is relatively robust and straightforward [194]. In addition to the model by Gosman et al. [192], there are other more complex models such as the Hall [195] or Legg & Raupach [196] model, which are particularly suitable for use in turbulent multiphase systems [194].
Bubble Aspect Ratio By default, a constant bubble diameter is modeled with a constant aspect ratio. To model a gas bubble size distribution, PBM are often used (Sect. 5.2). The larger the gas bubbles in the system, the more critical it is to model the aspect ratio. Different models have been developed. Besagni et al. [197] were able to show that the model by Wellek et al. [198] overestimates the aspect ratio compared to experimental data. The model of Vakhrushev & Efremov [199], which calculates the aspect ratio as a function of the Taylor number, showed good agreement with experimental data [199].

Forced-Aerated Bioreactors
The Euler-Euler model is the most widely used two-phase model for the modeling of forced-aerated bioreactors [107]. An overview of Euler-Euler simulations for bioreactors is given in Tab. 4, where all authors have calculated the k L a value. It can be concluded that the use of PBM is already standard and that Ansys Fluent is a popular tool for CFD modeling. An extensive overview of aerated and stirred bioreactors can be found in [176]. In general, a good agreement with experimental data can be achieved. For example, Bannari et al. [200] were able to show, depending on the process parameters, that the k L a values from simulation and experiment only differ by 6.8 % (without mixing) and 7.5 % (with mixing). The k L a values in the simulations of Azargoshasb et al. [201] differed from the experiments by 13 % at 600 rpm, 12.5 % at 800 rpm and 11 % at 1000 rpm. Amer et al. [17] were able to show in a 50-L single-use STR bioreactor that the use of PBM is necessary to perform a realistic simulation. The k L a value without PBM was 18.2 h -1 and with PBM 8.0 h -1 , compared to the measured value of 7.3 h -1 . However, Kaiser et al. [202] were able to achieve a deviation of only 2 % in the benchtop format without PBM at low aeration rates. At higher aeration rates, however, simulation and experiments deviated by up to 30 %.

Euler-Lagrange Approach
In contrast to the Euler-Euler approach, the Euler-Lagrange approach describes the disperse phase statistically. The individual particle paths are described by ordinary differential equations [220]. The use of the Euler-Lagrange approach for the two-phase modeling of STR is quite rare [221]. However, Wutz et al. [222] were able to show under different process conditions that k L a values for bioreactors in the field of mammalian cell culture can be successfully calculated using this approach. Weber and Bart [223] also used this approach for the modeling of bubble columns and were able to show that the results of the Euler-Lagrange and Euler-Euler approach are comparable. The main disadvantage of this approach is that the simulation time increases linearly with the number of bubbles [223]. This fact makes the Euler-Lagrange approach with the current computing power uninteresting for the economic characterization of bioreactors with high aeration rates as used, e.g., for microbial cultivation.

Lattice Boltzmann Approach
A relatively new method for two-phase CFD simulations is the Lattice-Boltzmann method (LBM) which, in contrast to the other methods presented here, does not numerically solve the Navier-Stokes equation but the Boltzmann equation. The method consists of two repetitive steps: the collision step and the streaming step. Further fundamentals can be found in the literature [224,225]. The use of LBM is currently still rare but seems to be a promising method for two-phase simulations [226]. The huge advantage of the LBM is the possibility of efficient parallelization [227]. Both Eibl et al. [228] and Sungkorn et al. [221] have used the LBM to model aerated STR. Shu et al. [229] were able successfully model a bubble column and validate the results, but no k L a value was calculated.

Population Balance Model
Population balance models (PBM) based on the population balance equation (PBE) can be used to describe the evolution of a population of particles [230]. One of the most widespread applications is the modeling of bubble size   distributions in dispersed gas-liquid systems [231]. For some years now, PBMs have been increasingly coupled with CFD models [232]. PBE describe the change of the number density function of one or more internal coordinates. In the case of gas bubble size distribution, typically the diameter of the gas bubble is considered, and the PBE corresponds to the Eq. (22) [232,233]. The source term consists of four different terms: death by coalescence D coal , birth by coalescence B coal , death by breakup D br , birth by breakup B br (Eqs. (23) to (26)). For the closure of the four source terms, coalescence and breakup models are used [113]. A broad overview of different breakup and coalescence models can be found in [96,101,234].
To simulate the bubble breakup, the breakup frequency is modeled on the one hand and the daughter size distribution resulting from the bubble breakup on the other hand, whereby most models only consider a binary bubble breakup [96]. Wang et al. [235] were able to show with experimental investigations that M-shaped daughter bubble size distributions as with the model of Lehr et al. [100] or Wang et al. [236] can achieve consistent results. The breakup model of Laakkonen et al. [237] was developed for stirred reactors and is very common in the literature. Other frequently used approaches are the models of Prince & Blanch [238], Luo & Svendsen [239] and Lehr et al. [100,139].
Coalescence is a more complicated phenomenon, since not only bubble-liquid interactions but also bubble-bubble interactions must be considered [97]. The coalescence frequency is modeled either by empirical models [240][241][242][243] or by the product of coalescence efficiency and collision frequency [100,235,238,[244][245][246]. Kaiser [202] was able to show in his investigations that the model by Lehr et al. has the lowest coalescence rate and, thus, overestimates the k L a value. Nevertheless, the model of Lehr et al. is one of the most widely used models [247]. Other frequently used models are those of Coulaloglou & Tavlarides [248], Prince & Blanch [238], Luo [249] and Wang [139,236].
To solve the PBE, a large number of solution methods have been implemented. An overview of conventional methods can be found in [250,251]. In commercial and open-source CFD codes, the fixed pivot method as a variant of the class method (CM, OpenFOAM 8) as well as variants of the method of moments (simulated method of moments and quadrature-based moment of methods (QBMM)) in Ansys Fluent 2020 and multi bubble size group (MUSIG) method in Ansys CFX 2020 are standard [252]. OpenQBMM, which is an extension of OpenFOAM, enables the use of the QBMM method for OpenFOAM, whereby good agreement with experiments has already been achieved in stirred aerated reactors [253]. In the class method, the gas bubble diameters of the gas bubble size distribution are divided into discrete classes. In order to obtain a solution that is not dependent on the choice of the number of classes, a sufficiently large number of classes must be chosen. Sanyal et al. [254] recommend at least 12-18 classes; other authors use 25 classes [255]. The most significant disadvantage of the class method is that the required computing time increases exponentially with the number of classes [254,256]. Another disadvantage is that gas bubbles cannot become larger than the defined maximum class diameter. However, due to the increasing calculation effort, a gas bubble maximum size must be defined [257].
Nevertheless, the maximum bubble size can be estimated according to Mersmann [258] for low viscous fluids. With methods of moment and its further developments, the gas bubble size distribution is not modeled as it is with the class method, but only the moments of the bubble size distribution. This leads to significantly shorter simulation times with the same accuracy [254]. The fact that only the moments of the gas bubble size distribution and not the distribution itself are known can be seen as the greatest disadvantage, but in practice often only the moments of distribution are of interest. For example, the Sauter diameter can be determined by the ratio of the third to the second moment [259].

Efficiency and Economy
Due to the nature of simulations, it is never possible to reproduce the exact reality with CFD simulations. CFD simulation should only be as accurate as necessary to be an economic complement to physical experiments. The errors that occur in a CFD simulation can be divided into model errors, discretization errors, iteration-convergence errors, rounding errors, as well as programming and user errors, whereby the latter are the only avoidable errors [18]. Model errors usually account for the largest fraction of the total error and are also difficult to estimate [18]. Concerning two-phase bioreactor simulations, the choice of a suitable two-phase model is of crucial importance (Fig. 2). An inappropriate model choice can lead to unrealistic results or unnecessarily increase the simulation time. For surface aerated systems, the use of a VOF model is therefore suitable. If forced-aerated bioreactors in benchtop scale with aeration rates as they are common in cell culture technology are to be modeled, the standard Euler-Euler model is entirely appropriate, since it could be shown that the influence of bubble coalescence and breakup has a negligible influence on the k L a value [202]. If the Euler-Euler model is used, additional models have to be taken into account, which also influence the model error and the simulation duration. These models include the drag force coefficient, lift force coefficient, virtual mass force coefficient, wall lubrication force coefficient, turbulent dispersion force coefficient, bubble aspect ratio and phase transfer. If, however, systems for microbial use or large bioreactors are investigated, CFD must be coupled with PBM for realistic modeling. In addition to the models mentioned previously, coalescence and breakup models are also used. For all three two-phase models mentioned here, a turbulence model must be chosen as well as models for fluid viscosity and surface tension.
A common method to estimate the discretization error is the Richardson extrapolation [260][261][262][263]. Other methods are the grid systematic refinement, grid convergence index, and the curve fitting method [263]. Apart from the discretization schemes, the mesh quality is of particular importance. A higher grid quality increases the accuracy of the approximation for surface and volume integrals and, thus, reduces the discretization error [264]. In order to carry out an economic simulation, the calculation grid should only be selected, so refined that the desired discretization error is achieved. With increasing grid density, the calculation time increases exponentially. However, with increasing grid density, the solution approaches the exact solution only asymptotically [14]. If the class method is used for solving the PBM, the computational time also increases with increasing class number (Sect. 5.2, population balance model).
Since the simulations are often transient, the time discretization must also be considered. To ensure stability of the transient simulation, the Courant-Friedrichs-Lewy number (CFL) is often used [265]. The CFL number must be less than or equal to one for explicit procedures to ensure stability. Implicit procedures can also be stable with CFL numbers greater than one, but the accuracy decreases with increasing CFL number. The CFL number fluctuates strongly over the whole grid if the grid is not adapted to the velocity. In this case, the maximum CFL number should be used to calculate the time step size [266]. Choosing an appropriate CFL number is challenging because a high CFL number is required to achieve rapid convergence. However, a low CFL number is required for reliable convergence [267].
To estimate the total error, it is necessary to validate the CFD simulations with physical experiments [268]. The measurement methods presented in Sect. 4.2 are suitable for validation with respect to the k L a value. If the flow field is also subject of validation, methods like particle image velocimetry [269][270][271][272][273][274], laser Doppler anemometry [275][276][277] or laser induced fluorescence [278,279] can be used. If an own CFD code is developed, a verification of the code must also be carried out beforehand [268].

Conclusions
Optimal oxygen supply for the cells in biopharmaceutical production processes is crucial, since both growth and product formation rate are directly dependent on it. It is therefore of interest to know the oxygen demand of the cells used and to determine the oxygen transfer rate of the system. With respect to the first point, it was possible to show that the specific oxygen uptake rates differ from organism to organism and the selected measuring principle. Empirical formulas can be used to estimate the k L a value. This is a fast method, but the values calculated can differ by a power of ten, depending on the model [256]. Furthermore, they do not offer any possibility of determining local limitations in the system. If a bioreactor is already physically existent, the k L a value can be determined experimentally by using the measuring methods described in Sect. 4.1. However, this can be very resource-intensive for large systems. CFD offers an excellent alternative to help solve this problem. It can bypass the inaccuracy of empirical correlations, is more resource-efficient when used correctly and offers the possibility of characterizing prototypes in silico while showing local limitations.
For surface-aerated bioreactors, the VOF approach provides a sufficiently good accuracy. Since gas dispersion rarely occurs in these systems, the consideration of forces between the two phases is negligible, and phenomena such as coalescence and bubble decay need not be taken into account.
For STR and airlift reactors, the VOF approach is unsatisfactory due to the oversimplification. The Euler-Euler approach is much better suited for disperse systems because the interfacial forces can be considered with this approach. To consider polydispersity and phenomena like coalescence and bubble breakup, the approach has to be coupled with PBM. This leads to more realistic results but extends the simulation time by a factor of approximately two, depending on the model. The Euler-Lagrange approach, which is less frequently used for the simulation of gas bubbles, is also suitable for such bioreactor systems. In the literature, equally good results are obtained, but the computational effort increases linearly with an increasing number of bubbles. In addition, with the Euler-Euler approach, a small increase can be observed in the computing time.
In the future, CFD will be used even more for the process engineering characterization of bioreactors. Due to the continuously increasing computing power with correspondingly decreasing prices, it will be possible to model phenomena like turbulence with more sophisticated models in economic timeframes. With available computing power, it is also possible to investigate a large number of process parameter combinations simultaneously. In addition, the coupling of CFD with growth and product formation kinetics will be increasingly used.
Stefan Seidel holds a bachelor's degree in biotechnology and a master's degree in applied computational life sciences from Zurich University of Applied Sciences. His focus is on the classical process engineering characterization of bioreactors as well as on multiphase CFD simulations and their validation using shadowgraphy and particle image velocimetry. Since 2017, he has been a research assistant in the Competence Center for Biochemical Engineering and Cell Cultivation Techniques at ZHAW and since 2020, a PhD student at the Technical University of Berlin. Valentin Jossen holds a PhD in process engineering from the Technical University of Berlin. Currently, Dr. Jossen is working at the Zurich University of Applied Sciences as a research associate in the Competence Center for Biochemical Engineering and Cell Cultivation Techniques. His research focus is on the characterization and development of new bioreactor systems, especially for the production of cell therapeutics, using classical process engineering and numerical approaches.
Dieter Eibl has held an engineering degree in food technology since 1981 and a PhD in biotechnology from the Technical University in Köthen since 1986. Since 1991, Prof. Eibl has been working at the University of Applied Sciences in Wädenswil as a lecturer. He is the head of the Competence Center for Biochemical Engineering and Cell Cultivation Techniques and working group leader for biochemical engineering. Prof. Eibl brings more than 25 years of professional expertise in upstream process development, scaling-up, project and team management.