**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*-*ε*-model (Eqs. (4) and (5)) and *k*-*ω*-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(\rho k)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_i}(\rho u\_i k) = \frac{\partial}{\partial \mathbf{x}\_i} + \left(\mu + \frac{\mu\_t}{\sigma\_k}\right) \frac{\partial k}{\partial \mathbf{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 \mathbf{x}\_i} - \rho \varepsilon \tag{4}$$

$$\frac{\partial(\rho\varepsilon)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_i}(\rho u\_i \varepsilon) = \frac{\partial}{\partial \mathbf{x}\_i} + \left(\mu + \frac{\mu\_t}{\sigma\_\varepsilon}\right) \frac{\partial \varepsilon}{\partial \mathbf{x}\_i} + \mathbf{C}\_1 \frac{\varepsilon}{k} \mu\_t \left(\frac{\partial u\_i}{\partial u\_j} + \frac{\partial u\_j}{\partial u\_i}\right) \frac{\partial u\_j}{\partial \mathbf{x}\_i} + \mathbf{C}\_2 \rho \frac{\varepsilon^2}{k} \tag{5}$$

*Computational Fluid Dynamics for Advanced Characterisation of Bioreactors Used… DOI: http://dx.doi.org/10.5772/intechopen.109848*

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 *Fp* and *F<sup>τ</sup>* 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)).

$$\mathcal{M} = \sum\_{i} \left( r\_i \times F\_{p\_i} + r\_i \times F\_{\tau\_i} \right) \tag{6}$$

$$F\_{p\_i} = p\_i \cdot A\_i \tag{7}$$

$$F\_{\tau\_i} = \mu \cdot \frac{\partial u\_i}{\partial \mathbf{x}\_i} \cdot A\_i \tag{8}$$

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 *ε* can be determined. This is either calculated directly, for example, when using the *k*-*ε*-model, or can be determined from *k* and *ω* when using the *k*-*ω*-model (eq. (9)). The energy dissipation rate can be used to predict the Kolmogorov length *λ<sup>k</sup>* (eq. (10)). If this is determined for each mesh cell, a Kolmogorov length distribution in the bioreactor can be determined [44].

$$
\varepsilon = \mathbb{k} \cdot w \cdot \boldsymbol{\beta}^\* \tag{9}
$$

$$
\lambda\_k = \left(\frac{\nu^3}{\varepsilon}\right)^{\frac{1}{4}} \tag{10}
$$

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*<sup>m</sup> can be fed into the system via a scalar transport equation (eq. (11)) [48].

$$\frac{\partial T\_{\rm m}}{\partial t} + \nabla \cdot (\mathbf{V} \cdot \mathbf{T} \mathbf{m}) - \nabla^2 ((D + D\_t) \cdot T\_{\rm m}) = \mathbf{0} \tag{11}$$

The subsequent transient observation shows how the tracer is distributed in the system by convection and diffusion (including turbulent diffusion *D*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 *χ* such as density *ρ* and viscosity *μ* are weighted by the phase fraction *α* (eqs. (14) and (15)). The surface tension *σ* can be considered here by using, for example, the continuum surface force (CSF) model [49].

$$\nabla \cdot \overrightarrow{v} = \mathbf{0} \tag{12}$$

$$\frac{\partial \rho \overrightarrow{\boldsymbol{v}}}{\partial t} + \nabla \cdot \left(\rho \overrightarrow{\boldsymbol{v} \, \overrightarrow{\boldsymbol{v}}}\right) = -\nabla p + \rho \overrightarrow{\mathbf{g}} + \nabla \cdot \nu \left(\nabla \overrightarrow{\boldsymbol{v}} + \left(\nabla \overrightarrow{\boldsymbol{v}}\right)^{T}\right) + \sigma \kappa \nabla a\_{\text{liquid}} \tag{13}$$

$$\chi = \sum \chi\_i a\_i, \quad \chi \in [\rho, \nu] \tag{14}$$

$$\sum a\_i = \mathbf{1} \quad \forall a\_i, \{a\_i | \mathbf{0} \le a\_i \le \mathbf{1}\} \tag{15}$$

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*L*a* value can also be determined for these systems, which employ pure surface-aeration. In contrast to experimental investigations, *k*<sup>L</sup> 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 *α*water ¼ 0*:*5. The *k*<sup>L</sup> 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*<sup>32</sup> ¼^ *d*) [53].

*Computational Fluid Dynamics for Advanced Characterisation of Bioreactors Used… DOI: http://dx.doi.org/10.5772/intechopen.109848*

$$k\_{\rm L} = \frac{2}{\sqrt{\pi}} \epsilon \nu \left(\frac{D\_{O\_2}}{\nu}\right)^{\frac{1}{4}} \tag{16}$$

$$a = \frac{A\_{\text{g}}}{V} = \frac{6a}{d\_{32}}\tag{17}$$

Seidel and Eibl [8] have been able to show that in a 3 L stirred bioreactor, under the same simulation conditions, the *k*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*<sup>L</sup> 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*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*S1 criterion [58–60]. However, CFD investigations concerning the *N*S1 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 g L�<sup>1</sup> (1% particle fraction) microcarriers, as Jossen et al. [69] did, would require approximately 40 � 106 carrier L�<sup>1</sup> 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–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.
