Computational Fluid Dynamics for Advanced Characterisation of Bioreactors Used in the Biopharmaceutical Industry – Part I: Literature Review

Abstract

Computational fluid dynamics (CFD) is a widely used tool for investigating fluid flows in bioreactors. It has been used in the biopharmaceutical industry for years and has established itself as an important tool for process engineering characterisation. As a result, CFD simulations are increasingly being used to complement classical process engineering investigations in the laboratory with spatially and temporally resolved results, or even replace them when laboratory investigations are not possible. Parameters that can be determined include the specific power input, Kolmogorov length, hydrodynamic stress, mixing time, oxygen transfer rate, and for cultivations with microcarriers, the NS1 criterion. In the first part of this series, a literature review illustrates how these parameters can be determined using CFD and how they can be validated experimentally. In addition, an overview of the hardware and software typically used for bioreactor characterisation will also be provided, including process engineering parameter investigations from the literature. In the second part of this series, the authors’ research results will be used to show how the process engineering characterisation of mechanically driven bioreactors for the biopharmaceutical industry (stirred, orbitally shaken, and wave-mixed) can be determined and validated using CFD.

1. Introduction

Bioreactors are closed systems that allow the cultivation of microbial, animal, and human cells in a controlled, aseptic environment. The simplest bioreactors are not instrumentalised (e.g., shake flasks) and dominate in mL-scale process development. In contrast, bench-top (10 L working volume), pilot (≥50 L working volume), and industrial-scale (working volumes of cubic metres) bioreactors are typically equipped with sensors, which allow regulation and adjustment of the cultivation process. At this time, the primary goal of such cultivation is the production of inocula, and bio-, cell-, or gene therapeutics. The bioreactor type and its operation parameters influence the growth and production behaviour of the production organism used. For example, the often-used microorganism Eschericha coli has high oxygen demands [1, 2] and a need for short mixing times [3] due to its fast metabolism. In addition, microorganisms such as E. coli are extremely robust thanks to their cell walls and small size. Therefore, the bioreactor can be operated with a high specific power input. Biopharmaceutical production with E. coli has a long tradition and dates back to 1982 when the production of insulin was approved by the Food and Drug Administration in the USA [4].

Today, most biopharmaceuticals, such as monoclonal antibodies, therapeutic hormones and many vaccines, are produced using mammalian cells, with Chinese hamster ovary (CHO) cells being the most frequently used [5]. Furthermore, human embryonic kidney (HEK) cells, Madin-Darby canine kidney (MDCK) cells, murine myeloma cells (NS0), and baby hamster kidney (BHK) cells are used [6]. The advantage of mammalian cells is that they can produce more complex molecules and can also perform glycosylation, which makes them even more attractive to use [7]. Furthermore, mammalian cells generally do not require as much oxygen as E. coli [1, 8]. However, they are larger and lack a cell wall, which makes them much more sensitive to fluid dynamic stress, a factor that has to be taken into account when selecting the bioreactor system and its operation parameters [9, 10].

In addition to mammalian cells and microorganisms, plant cells, insect cells and, more recently, stem cells are also used in biotechnological manufacturing processes. These cells are either continuous cells (unlimited life span) or primary cells (limited life span) and grow adherently (they need a surface to attach to, such as planar plastic surface, membranes or microcarriers) or in suspension.

Bioreactors can be classified in different ways. If the reusability of the cultivation container is taken into account, a classification into single-use and reusable systems can be made [11, 12]. If the type of mass transfer is taken into account, they can be categorised into static and dynamic systems, the latter of which can be further divided into mechanically, pneumatically and hydraulically driven bioreactors, depending on the type of power input [13, 14]. Mechanically driven systems, which are either stirred, orbitally shaken, or wave-mixed, predominate [13], the last of which are exclusively available as single-use variants. Single-use bioreactors have become very popular over the past 20 years [15]. The cultivation containers (rigid plastic vessels or flexible bags) are pre-sterilised, can be used immediately, and are disposed of after the cultivation has been completed. This makes the cultivation safer, helps to save time, may reduce costs, and often contributes to a lower impact on the environment despite the plastic waste that is generated. Efforts by manufacturers of single-use bioreactors have helped make them more user-friendly, for example, by replacing weak points such as the film material of the first and second bioreactor generation with improved ones that are less problematic in terms of leachables and extractables [16]. Currently, single-use bioreactors up to 6 m³ working volume are available on the market [15].

Stirred bioreactors (Figure 1a) are characterised by a mostly cylindrical vessel containing a stirring shaft with one or more stirrers. The dimensioning of the vessel as well as the type, number and positioning of the stirrer(s), which introduce power into the system, differ depending on the application. The oxygen required by the cells is added via the fluid surface or by active gassing with a sparger, depending on the application. Stirred bioreactors range from the mL-scale in research and development to more than 100 m³ in production processes [12]. In orbitally shaken systems, such as the shake flasks (Figure 1b), power is introduced into the system through the walls of the vessel by moving the entire system orbitally on a platform, and oxygen is typically only introduced via the free surface. Wave-mixed systems (Figure 1c) are also surface-aerated however, in these systems power is introduced by tilting the cultivation bag to create one or two degrees of freedom (DOF) of motion, resulting in an undulation in the fluid inside the bioreactor.

Figure 1.

Computer-generated images (CGI) of common mechanically driven bioreactor systems. a) Stirred bioreactor where the power is brought into the system by rotation of the stirring shaft and the stirrer (rotation marked by green arrow). b) Shake flask, which is an orbitally shaken bioreactor. Here the power is brought into the system by the orbitally shaking of the platform (movement of the platform marked by green arrow). c) Wave-mixed bioreactor, in which the power is brought in by the back and forth tilting of the bioreactor (rotation around tilting axis marked by green arrow).

Bioreactors should be both, technically characterised and biologically qualified. Details on biological qualification of bioreactors are described by Schirmer et al. [17, 18]. Technical or so-called classical process engineering characterisation includes the determination of specific power input, mixing or residence time, and oxygen transfer [19]. If adherent cells such as mesenchymal and induced pluripotent stem cells are grown on microcarriers, investigation of the suspension behaviour can also be beneficial [20]. Characterisation is typically performed in the laboratory by means of experiments (Section 4). However, such characterisation is resource-intensive (financial and labour costs). Moreover, only bioreactor systems that already physically exist can be characterised. However, with the help of computational fluid dynamics (CFD), it is possible to eliminate some of these disadvantages, making it is possible to characterise and optimise digital systems prior to construction. This allows expensive production stops to be avoided, and more complex investigations to be performed. Nevertheless, there are certain disadvantages associated with the use of CFD. Well-trained personnel and a lot of computing power are required. Validation is necessary to obtain reliable results, since CFD is only a model of the real world that is limited by a finite degree of accuracy. In this chapter, process parameters and their determination by CFD are described and discussed. In addition, possible hardware and software solutions are also presented as well as validation methods.

2. Process engineering characterisation

The specific power input P/V is one of the central process engineering parameters. It describes how much power per volume is introduced into the bioreactor system. The power input influences various other parameters such as oxygen transfer, mixing intensity and fluid dynamic stress [21]. As already mentioned, a distinction can be made between mechanical, pneumatic, and hydraulic power inputs. Typical specific power inputs for microbial cultivations are >5 kW m⁻³ [22] and thus significantly higher than for mammalian cells (5 to 310 W m⁻³) [23, 24]. The specific power input, which is mechanically introduced into the bioreactor system can be determined by the torque M at a given stirrer speed n for volume V (eq. (1)).

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

The specific power input determined by this method corresponds to the average power input. The average power input is equal to the product of the average energy dissipation rate $\overline{\epsilon }$ and the density $\rho$ (eq. (2)) [9]. However, the energy dissipation rate in a bioreactor is not uniform and can differ by several orders of magnitude. For example, in a stirred bioreactor, the energy dissipation rate is maximum at the stirrer tips and decreases sharply as the distance from the tips increases. These differences in energy dissipation rate can be expressed in terms of hydrodynamic heterogeneity, which is the quotient of ${\epsilon }_{\max }$ and $\overline{\epsilon }$. Orbitally shaken systems are characterised by high homogeneity, where ${\epsilon }_{\max }$/$\overline{\epsilon }$ ranges from 1 to 18 [25, 26], and wave-mixed bioreactors lie between orbitally shaken and stirred bioreactors (8.8 to 32.0) [27]. The greatest heterogeneity can be observed in stirred bioreactors, with values of up to 66 for 6-blade pitched turbines and as high as 147 for Rushton turbines [28]. Kolmogorov’s microscale theory, which is often used as an indicator of possible cell damage, can be derived directly from the turbulent energy dissipation rate. Assuming isotropic turbulence, this theory states that eddies break down into smaller and smaller eddies and dissipate into heat at the Kolmogorov length ${\lambda }_k$. Eddies that are significantly larger than cells (e.g., $20\mathrm{\mu m}$ for CHO cells) damage the cells less and simply transport the cells with them. However, if the vortices are in the size range of the cells or smaller, this can lead to damage to the cell membrane. In addition to the Kolmogorov length, fluid dynamic stress can also be used to assess possible cell damage. This is composed of normal $\sigma$ and shear $\tau$ stresses. Lethal and sublethal hydrodynamic stress varies greatly from one organism to another and can even vary for the same cell type [29]. In addition to the cell type, culture media composition [30, 31] and exposure duration also play crucial roles.

$$P/V=\frac{\sum \left({\epsilon }_i\right)V_i\rho }{V}=\overline{\epsilon }\cdot \rho$$E2

The mixing time is another process engineering parameter which is of importance for large bioreactors. The mixing time corresponds to the time required to achieve a defined mixing quality $M$ (typically 95% in biotechnology, eq. (3)). The mixing time should be kept as low as possible so that no oxygen and nutrient limitations or pH fluctuations occur. However, few exact values are published in the literature. For example, Anane et al. [32] showed that a clear metabolic switch occurred in a CHO cell line when the mixing time was longer than 90 s. However, the maximum mixing time suggested by Löffelholz et al. [33] for litre-scale devices is 30 s, which is below this critical value. A detailed description of the decolourisation method, which can be used to determine the mixing time, can be found in Bauer et al. [19] or in Maschke et al. [10].

$$M\left(t\right)=1-\frac{\mid c\left(t\right)-c_{\infty }\mid }{c_{\infty }}$$E3

A sufficient supply of oxygen is of critical importance in aerobic cultivations. In contrast to other nutrients, oxygen must be continuously added to the system, since it does not dissolve well in water-like media (Newtonian culture broths). The amount of oxygen required for a cultivation can be determined by the oxygen uptake rate (OUR), which is the product of the cell density and the specific oxygen uptake rate $q_{O_2}$. An overview of typical specific oxygen uptake rates can be found in Maschke et al. [10] and Seidel et al. [2]. To ensure sufficient oxygen supply, the oxygen transfer rate (OTR) must be equal to or greater than the OUR. The OTR is the product of the volumetric oxygen mass transfer coefficient $k_{\mathrm{L}}a$ and the difference between the oxygen saturation concentration $C^{\ast }$ and the dissolved oxygen concentration in the bulk liquid $C$.

For process characterisation purposes, the $k_{\mathrm{L}}a$ value is of particular interest. This value is the product of the liquid side mass transfer coefficient $k_{\mathrm{L}}$ and the specific interface $a$. In practice, the $k_{\mathrm{L}}a$ value is typically determined directly. The most common methods are the gassing-out method, sulfite method and the respiratory gassing-out method, the latter of which is a biotic method. A detailed description of methods for determining the $k_{\mathrm{L}}a$ value can be found in Seidel et al. [2] or Bauer et al. [19]. However, there are also measurement methods that allow the determination of just the specific interface (Section 4).

Another important parameter is the suspension criterion $N_{\mathrm{S}1}$. The characterisation of the sedimentation behaviour is particularly important for the cultivation of cells growing adherently on microcarriers, such as stem cells, which are also sensitive to hydrodynamic stress. Various authors have used the $N_{\mathrm{S}1}$ criterion to successfully grow stem cells on microcarriers [34, 35, 36]. This criterion describes the stirrer speed and thus the minimum power input required to completely suspend all microcarriers. It should be noted that at this speed only a suspension of the microcarriers is guaranteed, but not a homogeneous distribution in the cultivation vessel. This criterion ensures that not only are the cells subjected to the minimum necessary hydrodynamic stress, but also that there is no nutrient limitation due to heap formation. The $N_{\mathrm{S}1}$ criterion can be determined using a camera without the need for extensive apparatus.

3. CFD models

CFD is used in various fields, including biotechnology, for research, development and education [37]. The most widely used method is the finite volume method that discretises the Navier–Stokes equations. In addition, this method, which is used in most commercial and open source CFD codes (e.g., Ansys Fluent and CFX, Siemens Simcenter Star-CCM+, OpenFOAM), the finite element method (Autodesk CFD) or the Lattice-Boltzmann method (M-Star CFD, OpenLB) are also used. In the latter, the Boltzmann equation is discretised instead of the Navier–Stokes equation [38]. In this chapter, only the finite volume method will be discussed. To describe a single-phase laminar flow in a bioreactor, the conservation equations for mass and momentum are discretised and solved iteratively. Single-phase considerations are only suitable for stirred bioreactors where a constant fluid surface can be assumed (no vortex formation).

For most biotechnological applications, a turbulent flow needs to be generated to ensure sufficient mixing for the cells to be suspended, for the oxygen to be dispersed and for the nutrients to be evenly distributed. This is, however, not the case for bioreactors used for cultivating cells that may be sensitive to hydrodynamic stress. Several methods exist to account for turbulent flow, with the Reynolds-averaged Navier–Stokes (RANS) approach being the most widely used due to its economic advantages. This approach fully describes the turbulence using an additional model and only calculates the average flux. There are a number of different RANS turbulence models. However eddy-viscosity-based models, such as the $k$-$\epsilon$-model (Eqs. (4) and (5)) and $k$-$\omega$-model (equations can be found in Wilcox [39]), are the most widely used [40, 41]. If the phenomenon of turbulence needs to be considered in more detail, the unsteady-RANS, Hybrid RANS large eddy simulation (LES), LES or direct numerical simulation (DNS) can be used with increasing computational complexity. In fact, DNS is so computationally intensive that it is only used for simple geometries in research [42]. A detailed overview of turbulence models is described in Wilcox [39] and Rodriguez [43].

$$\frac{\partial \left(\mathit{\rho k}\right)}{\partial t}+\frac{\partial }{\partial x_i}\left(\rho u_{i}k\right)=\frac{\partial }{\partial x_i}+\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}+{\mu }_t\left(\frac{\partial u_i}{\partial u_j}+\frac{\partial u_j}{\partial u_i}\right)\frac{\partial u_j}{\partial x_i}-\mathit{\rho \epsilon }$$E4

$$\frac{\partial \left(\mathit{\rho \epsilon }\right)}{\partial t}+\frac{\partial }{\partial x_i}\left(\rho u_i\epsilon \right)=\frac{\partial }{\partial x_i}+\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_i}+C_1\frac{\epsilon }{k}{\mu }_t\left(\frac{\partial u_i}{\partial u_j}+\frac{\partial u_j}{\partial u_i}\right)\frac{\partial u_j}{\partial x_i}+C_2\rho \frac{{\epsilon }^2}{k}$$E5

Mechanically driven bioreactors are characterised according to the method used to introduce power into the system, described in Figure 1. For stirred systems, two different methods are mainly used to model the motion of the stirrer. With the frozen rotor method (also known as multiple reference frame method) the stirrer is not rotated, but in the zone around the stirrer the conservation equations are solved in a rotating reference system. Thus, there is no displacement or recalculation of the computational mesh, which makes this method resource efficient. However, this method is only suitable for steady-state considerations of stirred bioreactors. Another method is the sliding mesh approach, in which the mesh of the stirrer and stirred zone is shifted relative to the rest of the bioreactor. This method is significantly more computationally intensive than the frozen rotor method, but allows for a transient approach. There are also two different methods for orbitally shaken and wave-mixed bioreactors, both of which use a transient approach. The computational mesh can be manipulated by rotation and/or translation to mimic motion [44]. A less common method is to manipulate only the direction of the forces acting within the bioreactor. For example, Zhu et al. [45] successfully modelled orbitally shaken systems using manipulation of centrifugal forces.

In terms of the process engineering characterisation, a single-phase simulation can be used to determine the flow field, as well as the non-aerated specific power input, with the power input being determined analogously to eq. (1). Instead of measuring the torque, it can be determined using eq. (6) [46], where $F_p$ and $F_{\tau }$ correspond to the pressure and viscous force, and $r$ to the distance from the mesh cell centre to the axis of rotation (Eqs. (7) and (8)).

$$M=\sum _i\left(r_i\times F_{p_i}+r_i\times F_{{\tau }_i}\right)$$E6

$$F_{p_i}=p_i\cdot A_i$$E7

$$F_{{\tau }_i}=\mu \cdot \frac{\partial u_i}{\partial x_i}\cdot A_i$$E8

Normal and shear gradients can be determined from the velocity field or its gradients. A detailed derivation is described in Wollny [47]. By using RANS turbulence models, the turbulent energy dissipation rate $\epsilon$ can be determined. This is either calculated directly, for example, when using the $k$-$\epsilon$-model, or can be determined from $k$ and $\omega$ when using the $k$-$\omega$-model (eq. (9)). The energy dissipation rate can be used to predict the Kolmogorov length ${\lambda }_k$ (eq. (10)). If this is determined for each mesh cell, a Kolmogorov length distribution in the bioreactor can be determined [44].

$$\epsilon =k\cdot \omega \cdot {\beta }^{\ast }$$E9

$${\lambda }_k={\left(\frac{{\nu }^3}{\epsilon }\right)}^{\frac{1}{4}}$$E10

The mixing time in stirred bioreactors can also be determined based on a stationary single-phase simulation. For this purpose, the converged stationary solution of the velocity field is used as a basis. Subsequently, a virtual tracer $T_{\mathrm{m}}$ can be fed into the system via a scalar transport equation (eq. (11)) [48].

$$\frac{\partial T_{\mathrm{m}}}{\partial t}+\nabla \cdot \left(\mathit{V\cdot },{\mathrm{T}\mathrm{m}}_{}\right)-{\nabla }^2\left(\left(D+D_{\mathrm{t}}\right)\cdot T_{\mathrm{m}}\right)=0$$E11

The subsequent transient observation shows how the tracer is distributed in the system by convection and diffusion (including turbulent diffusion $D_{\mathrm{t}}$) over time. Analogously to the decolourisation method, the mixing time can thus be determined.

If the fluid flow in orbitally shaken, wave-mixed or vortex-forming stirred bioreactor systems is being investigated, the gas phase has to be modelled in addition to the liquid phase. Different methods exist for the investigation of muliphase systems. However, the volume of fluid (VOF) method is suitable for continuous two-phase systems, in which the continuity can be assumed. The continuity equation corresponds to eq. (12) and the momentum equation to eq. (13). The VOF method uses a mixed fluid in which the physical properties $\chi$ such as density $\rho$ and viscosity $\mu$ are weighted by the phase fraction $\alpha$ (eqs. (14) and (15)). The surface tension $\sigma$ can be considered here by using, for example, the continuum surface force (CSF) model [49].

$$\nabla \cdot \overrightarrow{v}=0$$E12

$$\frac{\partial \rho \overrightarrow{v}}{\partial t}+\nabla \cdot \left(\rho \overrightarrow{v}\overrightarrow{v}\right)=-\nabla p+\rho \overrightarrow{g}+\nabla \cdot \nu \left(\nabla \overrightarrow{v}+{\left(\nabla \overrightarrow{v}\right)}^T\right)+\mathit{\sigma \kappa }\nabla {\alpha }_{\text{liquid}}$$E13

$$\chi =\sum {\chi }_i{\alpha }_i,\chi \in \left[\rho ,\nu \right]$$E14

$$\sum {\alpha }_i=1\forall {\alpha }_i,\left\{{\alpha }_i|0\le {\alpha }_i\le 1\right\}$$E15

Similarly to the single-phase simulation, the power input can also be determined using the torque value. Here, the phase fraction must be taken into account as well as the fact that this is a transient determination, which may have to be averaged over time for the average specific power input. Shear and normal gradients as well as Kolmogorov lengths can also all be determined in a similar way to the single-phase simulations.

The $k_{\mathrm{L}}a$ value can also be determined for these systems, which employ pure surface-aeration. In contrast to experimental investigations, $k_{\mathrm{L}}$ and $a$ are determined individually in CFD simulations. For surface-aerated systems, the specific interface can be determined from the quotient of the liquid surface area and the working volume, with the liquid surface area determined using an iso-surface with ${\alpha }_{\text{water}}=0.5$. The $k_{\mathrm{L}}$ value can be approximated by several models. Higbie’s [50] model is widely used (eq. (16)) and requires the energy dissipation rate to be determined. The mixing time can be determined by a virtual tracer. However, in VOF models it is important to ensure that the virtual tracer does not migrate into the gas phase. There are different model approaches used in commercial and open source CFD codes (phasescalartransport in OpenFOAM, user-defined scalar in Ansys Fluent).

If, however, a dispersed system is to be considered, the VOF model is not able to represent this complex phenomena. Dispersed systems include active aerated systems and systems with particles such as microcarriers. In contrast to the VOF model, the Euler-Euler model does not use a mixed fluid, but instead solves the conservation equations for each phase individually. However, only one pressure field is calculated for all phases. Nevertheless, with the Euler-Euler model, interfacial forces such as drag, lift and virtual mass force, can also be taken into account, which makes it possible to model particles of a defined size (different particle sizes can be modelled by additional phases, but this is a computationally intensive method). These forces and their computations are described in detail in Seidel et al. [2], Pourtousi et al. [51] and Chuang and Hibiki [52]. If a known constant bubble diameter $d$ can be assumed, the specific interface $a$ can also be determined with the Euler-Euler model (eq. (17), for mono-disperse systems: $d_{32}\widehat{=}d$) [53].

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

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

Seidel and Eibl [8] have been able to show that in a 3 L stirred bioreactor, under the same simulation conditions, the $k_{\mathrm{L}}a$ value can be described as a power function of the bubble diameter with a negative exponent. However, combining CFD with population balance modelling (PBM) can minimise the heavy dependence of the simulated $k_{\mathrm{L}}$ value on the selected bubble diameter. This is achieved by modelling a population of bubbles with a bubble size distribution. In addition, PBM can also be used to model complex phenomena such as bubble breakup and coalescence by describing the change in the number density function of the gas bubble diameters, with the source term consisting of four components: death by coalescence, birth by coalescence, death by breakup, birth by breakup. A detailed description and overview of closure models for the population balance equations can be found in Liao and Lucas [54, 55], and a description of how the choice of coalescence and breakup models affects the $k_{\mathrm{L}}a$ value is provided by Seidel and Eibl [8]. Among the most frequently used approaches are the class method, which is used in OpenFOAM, and the method of moments, which is used in Ansys Fluent. A detailed overview of possible methods can be found in Nguyen et al. [56] and Li et al. [57], and their advantages and disadvantages are discussed in Seidel et al. [2].

The Euler-Euler model can also be used to perform suspension investigations and to determine the $N_{\mathrm{S}1}$ criterion [58, 59, 60]. However, CFD investigations concerning the $N_{\mathrm{S}1}$ criterion are not very widespread. In order to model the microcarriers as a solid phase, a Granular Flow Model (GFM) is required, with the kinetic theory of granular flow (KTGF) approach being widely used for this purpose [61, 62]. For modelling purposes, a granular pressure is introduced, which requires a gradient term for the granular phase momentum equation. In order to calculate the granular pressure, it is necessary to also calculate the granular temperature, which influences the collision behaviour of the particles [63]. An example of how the granular pressure can be modelled can be found in Syamlal et al. [64]. When using the Euler-Euler model, instead of modelling the individual particles (gas bubbles or microcarriers), only their phase fractions are modelled. If, however, the individual particles need to be modelled, the Euler-Lagrange model should be used, which describes the paths of the particles using an ordinary differential equation [65]. However, the Euler-Lagrange model does not consider the influence of turbulent flows on the particle path, which can be accounted for using the discrete random walk model [66]. It is a fact that the Euler-Lagrange model is used less frequently for biotechnological applications than the Euler-Euler model [2]. The main reason for this is that the computation time increases linearly as particle numbers increase [67]. Zieringer and Takors [68] recommend a maximum particle fraction of $10\%$ for the use of the Euler-Lagrange model. Accordingly, using $10\mathrm{g}{\mathrm{L}}^{-1}$ ($1\%$ particle fraction) microcarriers, as Jossen et al. [69] did, would require approximately $40\cdot {10}^6\text{carrier}{\mathrm{L}}^{-1}$ to be individually modelled [70].

Depending on the bioreactor system and the process parameters being investigated, different models should be considered in order to avoid uneconomical CFD simulations. It should be noted that irrespective of the model used for a CFD simulation, it will never be a true representation of reality. In addition to modelling errors, discretisation errors, iteration-convergence errors, rounding errors, programming errors and user errors can all lead to inaccuracies. Rounding errors are typically negligible (double-precision floating point numbers) and user errors should ideally be entirely avoidable. Discretisation error arise from the spatial discretisation of the mesh cells and the temporal discretisation of the time step sizes. A widely used approach for estimating discretisation errors arising from the choice of mesh is the Richardson extrapolation or the associated grid convergence index [71, 72, 73], which considers the change in the target parameter (e.g. specific power input or mixing time) using meshes of different resolutions. A detailed description of this procedure is described in Ramírez et al. [74] and Seidel et al. [44]. Alternatively, methods like the grid systematic refinement or the curve fitting method can also be used [75].

If the estimated discretisation error is known, an economically acceptable mesh can be chosen, since the computation time increases exponentially as the number of mesh cells increases, while the discretisation error approaches zero asymptotically [76]. For transient simulations, both spatial and temporal discretisation play a role. If the chosen time step size is too long, this can lead to instability and inaccuracies in the simulation. Like finer computational meshes, computing time also increases linearly as the time steps becomes shorter. The Courant-Friedrichs-Lewy (CFL) number helps to estimate the time step size [77]. If an explicit time integration scheme is used, CFL values above 1 can lead to non-physical values and divergence. A CFL number below 1 is preferable even in implicit schemes in order to avoid large errors. The CFL number, which is calculated individually for each mesh cell, can vary significantly within a single bioreactor model. Therefore, the maximum CFL number of any single mesh cell should not exceed 1 [78]. In order to determine the total error of the CFD results, validation by means of physical experiments is necessary. The most common validation procedures for CFD simulations of bioreactors are described in section 4.

4. Validation

The validation of CFD simulations is of great importance and is defined by the American Institute of Aeronautics and Astronautics [79] as “The process of determining the degree to which a model is an accurate representation of the real world from the perspective of the intended uses of the model.” In order to perform a validation, different experimental investigations can be carried out. It should be noted that the experiments are also not free of measurement errors, which must also be taken into account in the statistical statement. Depending on the target process variable, different methods are suitable for validating different bioreactor simulations. One of the most fundamental factors is the flow pattern. Particle image velocimetry (PIV) [80, 81, 82], laser doppler anemometry [83, 84] and laser induced fluorescence [81, 85] are suitable for measuring the velocities in a bioreactor and are especially used for stirred bioreactors. For orbitally shaken [86, 87] and wave-mixed [88] systems, often only the liquid height or distribution is measured.

If the power input is of relevance, it is reasonable to use it as a validation criterion. Thus, the measurement of the torque around the rotation axis is suitable for stirred [74, 89], orbitally shaken [90, 91] and wave-mixed [92] bioreactors. Alternatively, the specific power input can be determined either by measuring the electric current supplied to the motor [93, 94] or calorimetrically [95, 96]. Villiger et al. [97] described a method, which allows the maximum hydrodynamic stress to be determined using poly(methyl methacrylate) aggregates [98], which is also suitable for validation purposes. In addition, the mixing time can also be used as a validation parameter. Depending on whether the mixing time is evaluated globally or locally in the CFD simulation, the decolourisation method or the conductivity measurement according to DECHEMA e.V. Working Group for Single-Use Technology are both suitable for validation [19]. However, these are not the only valid approaches, for example, Vivek et al. [99] used a Raman spectroscopy probe to experimentally determine the mixing time to validate CFD-determined mixing times.

For aerated systems, either the $k_{\mathrm{L}}a$ value can be used directly as a validation parameter, or the bubble size distribution can be used. The $k_{\mathrm{L}}a$ value can be measured using different methods described in Seidel and Eibl [8] or Garcia-Ochoa and Gomez [100], and the bubble size distrbution by shadowgraphy [101], optical multimode online probe [102], capillary suction probe technique [103, 104], SOPAT endoscopy [105, 106] or focused beam reflectance [107, 108] measurements. The gas holdup, which is typically determined optically, by the differential pressure method [109] or radar level gauge [93] is also suitable for aerated systems.

For suspension studies with microcarriers or other particles, it is advisable to use the suspension behaviour or the $N_{\mathrm{S}1}$ criterion for validation. For this purpose, the particles that are to be suspended are added to an unstirred bioreactor. Then the speed is increased stepwise until the $N_{\mathrm{S}1}$ criterion is achieved. The evaluation is either performed by eye [110] or recorded with a camera [59, 111, 112, 113]. Delafosse et al. [113] and Loubière et al. [59] have refined this evaluation method by measuring the light attenuation using a light source mounted on one side of the bioreactor and the camera on the opposite side. Most authors have recorded lateral images as well as images from below, using a mirror for the latter. For the validation of their CFD simulations, Zhang et al. [114] recorded images with a high-speed camera. However, they did not examine cells grown on microcarriers, but rather immobilised Lactobacillus cells.

5. Hardware and software

In order to obtain a final result, several steps are necessary, which can be divided into preprocessing, processing and postprocessing [76, 115]. The individual steps are explained in more detail in Figure 2. Typically, several of these steps are included in today’s CFD software solutions. Using the Ansys software solution, the geometry can be created (SpaceClaim), the mesh can be generated (Meshing), the calculations can be performed (Fluent) and with CFD-Post the postprocessing can be completed. However, separate software is often used for geometry generation, such as Autodesk Inventor, SolidWorks, Solid Edge or Salome, and for postprocessing steps such as Tecplot and Paraview. The most widely used CFD software solutions for bioreactor modelling are Ansys Fluent and CFX (commercial) and OpenFOAM (open source). However, Simcenter Star-CCM+ from Siemens, Autodesk CFD, COMSOL Multiphysics and M-Star CFD are also used, with the latter specifically advertised for bioreactor applications (all commercial software).

Figure 2.

Visualisation of the process steps ranging from problem definition to validation. The Minifors 2 stirred bioreactor from Infors AG hase been used as an example.

In addition to the software, the hardware is also of critical importance for a simulation to be economic [116], with the hardware performance, the purchase price and the power consumption all playing a role. All current CFD software solutions allow for the parallelisation of the calculations. For this purpose, the computational mesh is partitioned into different domains. Individual processors then execute the computations in the individual domains. The communication between the domains is regulated, as usual in parallel computing, by means of the message passing interface (MPI) or other interfaces [116]. Different algorithms exist for creating partitions, which differ in their degree of automation and partitioning time [117]. For example, the Scotch algorithm of Pellegrini and Roman [118], which is based on dual recursive bi-partitioning and implemented in OpenFOAM, only needs the number of domains to perform the process. However, this algorithm requires a longer partitioning time than less automated algorithms [117].

While parallelisation reduces the required computing time, higher degrees of parallelisation also require more data to be exchanged between the processors, which has a negative impact on computing time. Figure 3 shows how the relative simulation time for modelling a stirred 3 L bioreactor changes as parallelisation increases (interFOAM, OpenFOAM v9). With the hardware setup described in Seidel and Eibl [8], the relative simulation time can be described as a power function with an exponent of $-0.817$. If there were no loss due to communication, an exponent of $-1$ would be expected. Harasek et al. [119] also performed parallelisation studies with up to 1024 cores, form which an exponent of $-0.93$ could be determined. According to Haddadi et al. [116], parallelisation should be performed in such a way that there are between 50,000 and 100,000 cells per domain, whereby the number of cells should be reduced as the complexity of the model increases.

Figure 3.

Relative simulation time depending on the processors used. A stirred bioreactor, which was examined with OpenFOAM v9, was used as a test case (VOF model, transient observation for 10 s). All simulations were performed three times to capture temporal variance. In each case, the Scotch algorithm was used for decomposition. The hardware setup is described in more detail in Seidel and Eibl [8].

Of the publications listed in Table 1, only one-third of the authors made statements about the hardware used. Between 4 [120] and 504 [121] cores were used for the calculations. Only five of the authors stated that they had used an HPC system [44, 121, 122, 123, 124], with the remaining authors having used desktop machines. The current versions of Ansys Fluent and Siemens Simcenter Star-CCM+ emphasised the use of graphics processing units (GPU) instead of central processing units for the calculations. Multiple GPUs can also be used here. Benchmarks from Ansys and Siemens show that using GPUs can shorten simulation time at the same time as reducing purchasing costs and power consumption [125, 126]. Of the authors listed in Table 1, only authors who used M-Star CFD stated that a GPU was used (M-Star CFD relies exclusively on GPUs). COMSOL, Autodesk CFD, Ansys CFX and OpenFOAM are not able to use GPUs for calculation by default. There are, however, a number of modifications such as MixIT, that allow OpenFOAM to us GPUs [127, 128]. Current developments show that classical CFD simulations could be complemented or replaced by other techniques. Machine learning techniques can be used to accelerate simulations or to improve turbulence modelling [129], and physics-informed neural networks (PINN) are increasingly used, for their ability to perform calculations 200 times faster to the same degree of accuracy [130, 131]. Another technique, which currently has no real application for bioreactor modelling, is quantum CFD (QCFD) [132, 133, 134, 135]. In quantum computing, the system can be in a superposition of multiple states at the same time. By adapting conventional algorithms to quantum computing, an enormous increase in speed could be achieved and the use of models for turbulence etc. would become obsolete [136].

Bioreactor [Configuration]	Parameter	Working volume	Software	Reference
Stirred tank bioreactors
Ambr15 and 250 (Sartorius AG) [2x3BSS, 2xRT6]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	15 mL to $250\mathrm{mL}$	Fluent, CFX	[120, 137]
Bioflo IIC (New Brunswick Scientific Co., Inc.) [2xRT6]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$4\mathrm{L}$	Fluent	[138]
Biotron GX (Biotron Inc.) [1xRT6, 1xST6, 1xP2]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$2.5\mathrm{L}$	Fluent	[139]
CellReady (Merck Millipore) [1xMI]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	1.5 L to $2.5\mathrm{L}$	Fluent	[140]
Commercially available single-use bioreactor [1xRT4]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$0.2{\mathrm{m}}^3$ to $1{\mathrm{m}}^3$	CFX	[141]
HyPerforma glass (Thermo Fisher Scientific Inc.) [1xPB4]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$2\mathrm{L}$	OpenFOAM	[8]
HyPerforma single-use bioreactor and Dynadrive (Thermo Fisher Scientific Inc.), n. S. [1x3BEE, 3xSD + 1xSD, 2x3BEE]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$2{\mathrm{m}}^3$ to $10{\mathrm{m}}^3$	Fluent	[121]
n. S. [2xRT]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$2\mathrm{L}$	Fluent	[142]
n. S. [1xPB6, 1xHPC, 1xRT6, 1xP2, 1xPB2]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$2\mathrm{L}$	Fluent	[143]
n. S. [1xRT6, 1xSD, 1xSD]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$10\mathrm{L}$	Fluent	[144]
n. S. [3xRT6]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$30\mathrm{L}$	CFX	[145]
n. S. [1xMI]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$50\mathrm{L}$	Fluent	[146]
n. S. [1xRT6 + 2x3BSS]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$0.2{\mathrm{m}}^3$	Fluent	[147]
n. S. (developed by Guoqiang Bioengineering Equipment Co.) [1xRT4]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	30 mL to 70 mL	CFX	[122]
BIOTECH-5JG (Shanghai Baoxing Bio-Engineer Equipment Co., Ltd.) [1xRT6]	$P/V$	$3\mathrm{L}$	Fluent	[148]
BioFlo 415 (Eppendorf SE) [2xRT6]	${\mathrm{\Theta }}_m$	$7\mathrm{L}$	CFX	[149]
CellReady (Merck Millipore), n. S. [1xMI]	${\mathrm{\Theta }}_m$	$2\mathrm{L}$ to $0.2{\mathrm{m}}^3$	Fluent	[150]
n. S. [1xMI]	${\mathrm{\Theta }}_m$	$50\mathrm{L}$	Fluent	[151]
n. S. [1xPB4, 1xRT2]	${\mathrm{\Theta }}_m$	$1.5\mathrm{L}$	OpenFOAM	[152]
BioFlo 110 (Eppendorf SE) [1xRT6 + 1x3BSS, 2xRT6]	${\mathrm{\Theta }}_m$, $P/V$	$2.2\mathrm{L}$	Autodesk CFD	[153]
n. S. [1xRT6, 1xPB6]	${\mathrm{\Theta }}_m$, $P/V$	$46.5\mathrm{L}$	Simcenter Star-CCM+	[154]
n. S. [1xRT8 + 1RT6]	${\mathrm{\Theta }}_m$, $P/V$	$54{\mathrm{m}}^3$	M-Star CFD	[155]
n. S. (Sartorius AG) [2x3BSS, 1xRT6 + 1x3BSS]	${\mathrm{\Theta }}_m$, $P/V$	$15\mathrm{L}$	Fluent	[123]
HyPerforma single-use bioreactor and Dynadrive (Thermo Fisher Scientific Inc.) [1xRT6 + 1x3BBS, 3xSD + 1xSD]	${\mathrm{\Theta }}_m$, $P/V$, $\tau$	1 L to 50 L	Fluent	[156]
BIOTECH-5JG (Shanghai Baoxing Bio-Engineer Equipment Co., Ltd.) [1xRT6]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, ${\mathrm{\Theta }}_m$	$5\mathrm{L}$	COMSOL; Multiphysics	[157]
n. S. [1x3BSS]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, ${\mathrm{\Theta }}_m$	$0.15{\mathrm{m}}^3$ to $0.35{\mathrm{m}}^3$	CFX	[158]
n. S. [4xRT6]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, ${\mathrm{\Theta }}_m$	$20\mathrm{L}$ to $100{\mathrm{m}}^3$	Fluent	[159]
XcellerexTM XDR-10 (Cytiva) [1xPB3]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, ${\mathrm{\Theta }}_m$	4.5 L to 10 L	OpenFOAM, MixIT	[160]
n. S. [4xRT6]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, ${\mathrm{\Theta }}_m$, $P/V$	$38\mathrm{L}$	Fluent	[161]
Xcellerex XDR-200, 500 and 2000 (Cytiva), Applikon 5 L (Applikon Biotechnology B.V.) [1xPB3]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, ${\mathrm{\Theta }}_m$, $P/V$	$3.9\mathrm{L}$ to $2{\mathrm{m}}^3$	M-Star CFD	[162]
n. S. [1xRT6+2x3BSS]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, $P/V$, ${\lambda }_k$	$0.3{\mathrm{m}}^3$	Fluent	[124]
n. S. [1xRT6, 2xRT6]	${\lambda }_k$	$3\mathrm{L}$ and $0.3{\mathrm{m}}^3$	Fluent	[163]
n. S. [1xA310, 1xRT6, 1xPB6]	Sedimentation	$24\mathrm{L}$	M-Star CFD	[164]
Wave-mixed bioreactors
Wave bioreactor 20/50 (GE Healthcare)	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$5\mathrm{L}$	Fluent	[88]
ReadyToProcess WAVE 25 bioreactor (Cytiva)	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, $P/V,\tau$	$2\mathrm{L}$ and $20\mathrm{L}$	OpenFOAM	[165]
CELL-tainer Utility (Celltainer Biotech BV) [with and w/o expansion channels]	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, ${\mathrm{\Theta }}_m$, $P/V$	0.15 L to 20 L	OpenFOAM	[44]
n. S. (Wave bag)	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, ${\mathrm{\Theta }}_m$, $P/V$	$10\mathrm{L}$	CFX	[27]
Cellbag CB2L and CB20L (Cytiva)	Sedimentation, $\tau$	$2\mathrm{L}$ and $20\mathrm{L}$	Fluent	[166]
Orbitally shaken bioreactors
TubeSpin 50 (Techno Plastic Products AG)	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	10 mL to 50 mL	Fluent, OpenFOAM	[167, 168]
Shake flask (a)	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	15 mL to 200 mL	Fluent, CFX, OpenFOAM	[87, 168, 169, 170]
TubeSpin 600 (Techno Plastic Products AG)	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	100 mL to 500 mL	Fluent	[45]
Orbitally shaken bag (Infors AG)	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	1 L to 10 L	OpenFOAM	[168]
n. S. (Cylinder)	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$	$10\mathrm{mL}$ to $0.1{\mathrm{m}}^3$	Fluent, OpenFOAM	[171, 172]
n. S. (Cylinder)	$P/V$	$10\mathrm{L}$ to $0.2{\mathrm{m}}^3$	OpenFOAM	[173]
TubeSpin 50 and 600 (Techno Plastic Products AG), Optimum Growth flask 500 and 5000 (Thomson Instrument Company), Shake flask	$P/V$	10 mL to 3000 mL	OpenFOAM	[10]
SB10-X (Adolf Kühner AG)	${\mathrm{k}}_{\mathrm{L}}\mathrm{a}$, $\tau$	2.5 L to 10 L	Fluent	[174]

Table 1.

Overview of process characterisations of bioreactors using CFD.

The bioreactors were subdivided according to their power input into stirred, wave-mixed and orbitally shaken systems. 3BEE, 3-Blade elephant ear stirrer; 3BSS, 3-Blade segment stirrer; A310, Hydrofoil A310; MI, Marine impeller; P, Paddle stirrer; PB, Pitched blade stirrer; RT, Rushton turbine; SD, Special design; ST, Smith turbine; the number after the label indicates how many stirrer blades were present.

6. Conclusions and outlook

CFD is a widely used tool that has played a role in biotechnology for many years and is constantly gaining in popularity. It can be used to determine classical parameters such as power input, mixing time and oxygen transfer, as well as to carry out further evaluations with spatial and time resolutions. CFD is also an excellent approach for determining hydrodynamic stress. However, depending on the bioreactor system and parameters, different modelling approaches should be selected to determine these parameters. Investigations have shown that for bioreactor modelling, the Fluent commercial software from Ansys dominates, followed by OpenFOAM, an open source solution. Very few authors compute on GPUs, however they still use CPUs. Even tough when used in conjunction with the latest software solutions, GPUs cost less, use less power and perform faster computations. In spite of the considerable computing power currently available, phenomena such as turbulence and gas bubble breakup in industry-relevant simulations, such as bioreactors, are still approximated using semi-empirical models. In the future, more powerful GPUs or even quantum processors for QCFD will increase computational power, and methods like PINN will speed up computations, potentially making it possible to eliminate the need for these models and perform even more accurate and faster simulations. In the second part of the series, the methods described here for process engineering characterisation are demonstrated based on real-world examples, using both commercial and open-source software.

Acknowledgments

We would like to thank Darren Mace for English proofreading.

Conflict of interest

The authors declare no conflict of interest.