Validation of Novel Lattice Boltzmann Large Eddy Simulations (LB LES) for Equipment Characterization in Biopharma

Abstract

Detailed process and equipment knowledge is crucial for the successful production of biopharmaceuticals. An essential part is the characterization of equipment for which Computational Fluid Dynamics (CFD) is an important tool. While the steady, Reynolds-averaged Navier–Stokes (RANS) k − ε approach has been extensively reviewed in the literature and may be used for fast equipment characterization in terms of power number determination, transient schemes have to be further investigated and validated to gain more detailed insights into flow patterns because they are the method of choice for mixing time simulations. Due to the availability of commercial solvers, such as M-Star CFD, Lattice Boltzmann simulations have recently become popular in the industry, as they are easy to set up and require relatively low computing power. However, extensive validation studies for transient Lattice Boltzmann Large Eddy Simulations (LB LES) are still missing. In this study, transient LB LES were applied to simulate a 3 L bioreactor system. The results were compared to novel 4D particle tracking (4D PTV) experiments, which resolve the motion of thousands of passive tracer particles on their journey through the bioreactor. Steady simulations for the determination of the power number followed a structured workflow, including grid studies and rotating reference frame volume studies, resulting in high prediction accuracy with less than 11% deviation, compared to experimental data. Likewise, deviations for the transient simulations were less than 10% after computational demand was reduced as a result of prior grid studies. The time averaged flow fields from LB LES were in good accordance with the novel 4D PTV data. Moreover, 4D PTV data enabled the validation of transient flow structures by analyzing Lagrangian particle trajectories. This enables a more detailed determination of mixing times and mass transfer as well as local exposure times of local velocity and shear stress peaks. For the purpose of standardization of common industry CFD models, steady RANS simulations for the 3 L vessel were included in this study as well.

1. Introduction

In the pharmaceutical industry, CFD is no longer only an innovative tool that has to prove its value but has evolved to be part of the day-to-day business. In fact, due to increasing computational power and more sophisticated simulation software, CFD is typically integrated in scale-up models, process design, optimization and the development of scale-down models [1,2,3,4]. Thereby, its applicability is shown to address issues in up- and downstream development as well as for the production of biopharmaceuticals. Typical topics are risk mitigation during transfer, troubleshooting and reduction in lab experiments [5,6]. In line with the digitalization trend in the biopharmaceutical industry, other modeling approaches for upstream and downstream processing have also recently gained recognition [4,7,8,9]. The benefits are reduction in costs and improved equipment and process understanding, which can ultimately lead to a faster time to market.

Several CFD methods and turbulence models exist. Especially for the simulation of mixing tanks, different approaches were proposed [10,11,12,13]. In general, no standard procedure was specified on how a simulation should be set up and which steps are necessary to verify and validate applied numerical codes. However, the U.S. Food and Drug Administration provides recommendations to the industry on the formatting, organization, and content of reports for computational fluid dynamics [14]. Nevertheless, a validation of conducted simulation studies is inevitably dependent on the examined problem.

Stirred vessels are commonly used in biopharmaceutical production due to their superior combination of mixing performance, mass and heat transfer, which are especially important for aerated processes [15]. These vessels are usually characterized by the hydraulic power introduced into the system. The power consumption can be calculated from the overall torque derived from the pressure and tangential stress distribution on steady and moving walls. The resulting power numbers were shown to be predicted accurately by the Reynolds-averaged Navier–Stokes (RANS) k − $\epsilon$model, using the steady-state Multiple Reference Frame (MRF) approach to account for impeller rotation [16,17,18]. This computationally inexpensive method was proven to be a valuable tool for efficient equipment characterization to facilitate scale-up, as the simulation time was reported to be reduced by an order of magnitude, compared to other turbulence models [19,20]. Numerical grid independence for power number determination was already shown for coarser meshes; however, the obtained power numbers were also influenced by the size of the chosen volume within the rotating reference frame [21,22].

For detailed process and equipment knowledge, the resolution of transient structures within a stirred tank is essential. In particular, the numerical determination of mixing times relies on the correct depiction of time-resolved flow fields to capture the inherently turbulent, non-stationary character of the flow [23]. According to Haringa et al. [24,25], transient simulation schemes, such as the RANS sliding mesh approach or Large Eddy Simulations, outperformed MRF in single-phase stirred tanks. Therefore, LES is the preferred choice if facilities allow for finer meshes, additionally needing no definition of an MRF. The increased computational demand can be covered more efficiently when using parallel-computing platforms for which the Lattice Boltzmann (LB) method is ideally suited. Runtime performance of multiple orders of magnitude faster can be achieved on desktop graphic processing units (GPU) in contrast to central processing unit (CPU) clusters. In combination with LES, good performance was shown, simultaneously allowing for a coarser grid size, compared to direct numerical simulations (DNS) [26,27].

Advantages and limitations of the RANS k − ε model MRF have been extensively discussed in the literature. Prediction quality of obtained steady flow fields and power numbers are shown by Particle Image Velocimetry (PIV) or laser-doppler anemometry (LDA) and torque measurements [11,12,19,28]. LB methods have just recently become popular in the biopharmaceutical industry due to the availability of new commercial solvers [27,29,30]. Only a few articles have been published on LB LES for stirred tanks, lacking information on the quality of flow field prediction. As this approach offers the possibility to examine transient flow structures, an appropriate experimental validation method is necessary. For this reason, a new approach was chosen, which considers the three-dimensional Lagrangian flow structures measured by particle tracking methods [31]. Conventional approaches based on statistical or time-averaged data, such as PIV or LDA, do not provide such detailed data. Due to very high spatial as well as temporal resolutions of up to 150 µm/px and 800 Hz respectively, even more extensive validation datasets can be provided. This is due to the fact that not only can the three-dimensional velocity components be highly resolved, but this information is also available for almost every point inside the stirred tank reactor.

Standards for simulation procedures for stirred tanks in industrial applications are still lacking, making it hard to compare results from different publications or companies [4]. In this study, we present two standard approaches for bioreactor characterization within an industrial environment. The state-of-the-art RANS k − ε model MRF steady simulation was used for fast and easy determination of power numbers, whereas high resolution LB LES were used to display transient flow structures. We conducted a grid refinement study as a standard in both cases to verify numerical independence of the solution. Presenting our internal equipment characterization procedures, steady simulations were extended by a study varying the volume of the rotating reference. Furthermore, for using LB LES, initial validation of transient flow fields at small scale was deemed necessary. Therefore, the comparison to novel four-dimensional Particle Tracking Velocimetry (4D PTV) measurements was implemented. Based on the 4D PTV measurements, not only can the transient and time-averaged flow fields now be used for the validation of the numerical simulation, but, for the first time, also the Lagrangian particle trajectories.

2. Materials and Methods

2.1. Reactor Setup

A laboratory scale stirred tank reactor with a working volume of 2.8 L is used for numerical simulation and 4D PTV measurements as shown in Figure 1. The tank with a diameter of D = 0.13 m, a liquid filled height of H = 0.228 m is equipped with two Rushton impellers (d = 0.036 m). The off-bottom clearance is $C$ = 1.5 d and impeller spacing $\mathsf{\Delta }C$= 2.2 d. Single phase experiments and simulations are performed in water ($\rho$ = 998.2 kg/m³, $\eta$ = 1.0016 mPas, $\sigma$ = 0.0728 N/m) at 20 °C. A baffled (3 baffles 120°) and an unbaffled system are examined. Stirring speed is set to n = 250 and 350 rpm. For better optical access to the stirred vessel, an octahedral water basin is installed around the stirred vessel in the experiment. In addition, the baffles are made of acrylic glass to eliminate shadowing by the baffles. Further details about the experiments will be published separately.

2.2. Experimental Setup

The experimental validation of the transient numerical simulations of the 3 L stirred tank reactor (descripted in Section 2.1) is based on spatial (three dimensions) and temporal (plus 1 dimension) high resolution 4D PTV measurements. Based on the observed Lagrangian particle trajectories, not only the classical two-dimensional flow fields can be resolved but rather the three-dimensional flow structure in the stirred tank is determined, which provides much deeper information for the validation of numerical simulations. The particle tracking method used is the “Shake-The-Box” approach, which is currently the most efficient particle-tracking algorithm available [31]. The acquisition of the images, as well as the application of the 4D PTV algorithms, is done using the commercial software FlowFit from LaVision (Göttingen, Germany).

To determine the Lagrangian particle motion 150–180 µm size fluorescent polystyrene particles are used (Cospheric LCC, Santa Barbara, CA, USA). A particle per pixel density of 0.005 ppp is aimed at in the experiments, which corresponds to approx. 5000 particles in the reactor volume. For the excitation of the fluorescence signal, three high-power LEDs are used, which are synchronized to the camera signal by means of a programmable timing unit (PTU) from LaVision. The particle movements are recorded for a time interval of t = 3.75 s by four high-speed cameras. The spatial resolution for this setup is s = 150 µm/px and the temporal resolution in between two recorded images is dt = 0.0125 s. The reactor volume, as well as the perspective, are geometrically calibrated by means of a traversable calibration target inside the reactor. In addition, a volume self-calibration and subsequent calculation of the optical transfer function (OTF) is performed to improve the geometric calibration and particle detection by the algorithm [31].

For the measurement of the power number, a motor (ViscoPaktRheo X7, HiTec Zang, Germany) with integrated torque sensor is used. A detailed description of the methodology and the data evaluation is given in [32].

2.3. Numerical Simulations

The pressure-based solver from the commercial finite volume-based software ANSYS Fluent 19.0 (ANSYS Inc., Canonsburg, PA, USA), using the RANS approach, is employed to perform steady simulations. The incompressible continuity equation is integrated over a control volume according to the following:

(1)$\nabla U=0$

with $U$ as the three-dimensional velocity vector (Reynolds averaged). The equation of momentum conservation can be written as follows:

(2)$\rho \left(\frac{\partial U}{\partial t}+U \nabla U\right)=-\nabla p+g \rho +\nabla \left(\eta \nabla U\right)-\rho \overline{U_i^{\prime }{U}_j^{\prime }}$

with $\rho$ as density, $p$ as pressure, $g$ as gravitational acceleration, $\eta$ as dynamic viscosity and $\rho \overline{U_i^{\prime }{U}_j^{\prime }}$ as Reynolds stresses consisting of the fluctuation velocities $U_i^{\prime }$ and $U_j^{\prime }$, which are modeled by the eddy viscosity principle after Boussinesq [33], using the realizable $k$ − $\epsilon$ turbulence model.

The enhanced wall treatment is enabled, and curvature correction is applied. The multiple-reference frame (MRF) model is used to account for impeller rotation. All walls are set to no-slip boundary conditions and the top is treated as a free-surface. The second order upwind scheme is used for the spatial discretization of momentum, turbulent kinetic energy and dissipation rate. In order to guarantee numerical independence of the obtained results, a mesh study is performed using different number of polyhedral cells ranging from 100,000 up to 700,000 with a fixed rotating reference frame volume (see Figure S1). The chosen impeller diameter to reference frame diameter (D-ratio) and the blade height to height of the rotating reference frame (RRF) volume ratio (H-ratio) are 1.5. As flow characteristics are significantly influenced by the size of the inner volume of the RRF within the MRF approach, the volume of the rotating frame is varied around the top and bottom impeller, respectively. For each condition, the agitation speed is varied to display the evolution of the power number over the stirrer speed. Convergence of the simulations is assumed if the residuals for all equations$\approx$10⁻⁴; however, for higher rotating speeds, only $\approx$10⁻³ was reached. The power number is calculated according to the following:

(3)$Po=\frac{P}{\rho n^3 d^5 }$

with $\rho$ as density, $n$ as stirrer frequency, $d$ as impeller diameter and the power

(4)$P_t=2 \pi M n$

with the torque $M$and the stirrer frequency $n$.

Transient simulations are conducted with the commercial software M-Star CFD 2.10.57 (M-Star Simulations, LLC, Ellicott City, MD, USA), which uses the Boltzmann transport equation to model the transport carrier probability density [34,35]:

(5)$\frac{\partial f}{\partial t }+U \nabla f=\mathsf{\Omega }\left(f,f\right)$

with $f$ being the probability density function, $U$ the three-dimensional velocity vector and $\mathsf{\Omega }$ as collision operator. By linking $f$ to the fluid momentum and $\mathsf{\Omega }\left(f,f\right)$ to the fluid kinematic viscosity, macroscopic variables are recovered, and the Navier–Stokes equations are solved. Therefore, the molecular velocities are discretized into the stable velocity vector set D3Q19 [36].

Turbulence is modeled using Large Eddy Simulation (LES). The effects of the filtering are controlled by the Smagorinsky–Lilly subgrid closure model with a Smagorinsky coefficient of 0.1 [37,38]. The single-phase model is applied for the baffled system and the free-surface model involving a moving free surface is used for the unbaffled system. A free slip boundary condition is assumed at the free surface. Static walls are treated with the no-slip, bounce-back approach [38] and interactions between the fluid lattice and the impeller are handled by the immersed boundary method [39]. Time step size $\mathsf{\Delta }t$ is chosen to obtain a Courant number $Co$ of 0.05 to keep lattice density fluctuations around 1%, with tip speed $V_{tip}$ as reference velocity with the following:

(6)$Co=V_{tip}\frac{\mathsf{\Delta }t}{\mathsf{\Delta }x}$

A mesh independence study is performed using different lattice spacings from 1.7 million up to 379 million lattice points, resulting in a grid spacing from $\mathsf{\Delta }x$ = 1.3 mm to $\mathsf{\Delta }x$ = 0.22 mm. Therefore, the time averaged radial and axial velocity profiles as well as the evolution of the Power number are compared. Time averaged quantities are recorded from the developed flow field. The Power number of the power determined by torque $P_t$ and by the following integral predicted energy dissipation:

(7)$P_{\epsilon }={{\displaystyle \int }}^{}\epsilon \rho dV$

are compared.

To compare simulation results with 4D PTV data, virtual particles are injected as massless tracers that follow the fluid streamlines after the flow field is converged. Particle movements are recorded, and flow variables are averaged over a time interval of 3.75 s with a sampling rate of dt = 0.0125 s according to the experiment.

3. Results and Discussion

A laboratory-scale stirred tank reactor is characterized numerically by two different approaches. The steady approach with ANSYS Fluent is used to determine the Power numbers and to demonstrate a fast and inexpensive equipment characterization. Transient simulations are conducted with M-Star CFD to gain more detailed knowledge about the flow fields in stirred tanks. Transient and mean flow structures are validated by means of 4D PTV data.

3.1. Steady State Simulations

A mesh study for the baffled 3 L system is performed to prove numerical independence of the obtained results. Therefore, the agitation rate is varied for each grid size within the specified range, and $P_t$ as well as $P_{\epsilon }$are determined by averaging over the turbulent regime (Re > 10⁴). As shown in Figure 2A, a mesh of about 350,000 polyhedral cells is fine enough to determine power numbers $Po$ by torque, as the power number converges to a stable value of $\approx 8.1$. Consequently, the rotating frame volume study is performed with this mesh. Convergence for the power number $Po$ at a relatively coarse grid size is also shown in the literature [16]. The influence of the rotating reference volume size on the evolution of the power number is clearly visible in Figure 2B. If the chosen impeller diameter to the reference frame diameter (D-ratio) and the blade height to the height of the rotating reference frame (RRF) volume ratio (H-ratio) is below 1.5, the value for the power number is increasing sharply. Whereas an elevated D-ratio showed no impact until a ratio of 2, higher D-ratios resulted in a decrease in the power number. Increasing the H-ratio further up to a H-ratio of 8 led to a convergence of power number up to a mean value of$Po=8.4\pm 0.2$. A combined volume of the upper and lower impeller results in a power number $Po$= 7.2 $\pm$ 0.2 for increasing height of the combined RRF (data not shown). Deviations in the flow field prediction dependent on the rotating reference frame size have been reported in several studies [22,40,41], leading to the general conclusion to enlarge the dimension of the rotating part toward a region where the flow variables’ gradient is low to reduce numerical errors due to interpolation inaccuracies in transferring the results between both domains. Some authors suggest locating the interface between the rotating and stationary frame 0.5 impeller radius from the impeller tip, whereas others recommend a ratio of 1.5 [21,40]. For this study, a D-ratio between 1.5 and 2.0 is considered appropriate, as higher RRF diameters would place the interface of the RRF too close to the baffles, where too-high velocity gradients exist. Additionally, for the axial extension of the rotating zone, different propositions are made reaching from one to four impeller radius above and below the impeller [42]. Although in this study, the value for the combined RRF is closer to the experimental value of $Po$ = 7.5, marked as a dashed line in Figure 2, the converged power number $Po$ for the increased H-ratios of the two separate RRF is preferred, as huge instabilities within the prediction of turbulent parameters are noticed for the combined RRF. Considering that in typical cell culture bioreactors, many installations, such as probes or sampling ports, are installed that influence the flow structure, the size of the rotating domain has to be evaluated for each individual case. Especially since the influence of the rotating domain size on power number determination is as significant as the grid size, we recommend a RRF size study for reliable power number determination.

3.2. Transient Simulations

Transient profiles are examined for the baffled system at two different stirring frequencies. Additionally, velocity profiles of the system without baffles are compared. To prove mesh convergence, a grid refinement study is performed showing the velocity profiles and the evolution of power number. Using the respective grid, time averaged and transient flow structures are compared to 4D PTV data.

3.2.1. Grid Refinement Study

The time-averaged profiles of the normalized velocity magnitude for the baffled system and a stirrer frequency of $n$= 250 rpm are depicted in Figure 3. The velocity is azimuthally averaged over the radial position within the impeller discharge stream of the respective impeller. In this case, numerical independence of the velocity is reached at 27 million grid points, as further refinement resulted in the same velocity profile. This is also reflected in the averaged axial profiles of the normalized velocity magnitude, which are displayed in Figure S2. In Figures S3 and S4, the averaged radial and axial profiles for an impeller frequency of 350 rpm are displayed showing a mesh independence for identical grid size. A grid refinement study for the unbaffled system is not considered, as smaller velocity gradients are to be expected and therefore, a grid resolution of the baffled systems is considered sufficient. A similar grid resolution was previously found to be adequate [27].

The evolution of power number $Po$ determined by torque $P_t$ and by integral predicted energy dissipation $P_{\epsilon }$ is depicted in Figure 4A. Only minor deviations between these two values reaffirm the self-consistency of the solver. With increasing grid density, the power number stabilizes to a value of $Po=$ 7.5 $\pm$ 0.5, which fits the experimentally obtained value of $Po=$ 7.5 marked as a dashed line. Figure 4B shows the experimental and numerical power number $Po$as function of the Reynolds number (Re). Simulated power numbers with two Rushton turbines are in good agreement with the results by Zlokarnik [43], considering that the authors used only one Rushton turbine for their experiments. If regimes are divided into a laminar regime with $Po$ being proportionally dependent on Re, a transient and a turbulent regime where a further increase in Re shows no impact on $Po$, the simulation results show that the LB LES approach is capable of depicting this trend not only for the turbulent regime. Consequently, the LB LES approach is applicable to various flow regimes.

3.2.2. Validation of Numerical Simulations by 4D PTV Data

Simulated time-averaged flow fields and transient structures are compared to 4D PTV measurements. For the experimental and time-averaged velocity fields, an interpolation of the particle trajectories onto a Eulerian grid is created, using the FlowFit software from LaVision. For this purpose, a subvolume size of 42 voxels with an overlap of 75% is chosen. To improve the statistical accuracy, at least 10 tracks must have passed a voxel in order to calculate a velocity vector for the corresponding voxel. The mean velocity magnitude of the baffled system at n = 250 rpm is depicted in Figure 5A as a vector field, which represents a slice through the reactor. For better visualization, the three-dimensional velocity vectors are shown as a projection onto the sliced plane. Typical flow structures of radial pumping Rushton turbines, showing the characteristic upper and lower recirculation loops around the impeller, are visible for both the simulation and experiment. Highest velocities close to the impeller tip speed of 0.47 m/s are obtained within the impeller discharge stream. In general, the experimentally measured flow field is represented very well by the simulated velocity profile, locally matching the value and orientation of the velocity magnitude. Very good results are also obtained for other conditions tested, with maximal velocities being slightly underpredicted in the experiment for n = 350 rpm (see Figure S5A) and smaller deviations of simulated and experimental velocity magnitudes around the impeller shaft for the unbaffled system (see Figure S5B). The time-averaged velocity magnitude distribution is shown in Figure 5B for two different agitation frequencies. One third of the bioreactor is evaluated due to the rotational symmetry and the best optical access achieved in the front part of the vessel. At $n$ = 250 rpm, velocity magnitude of $V_{mag}$ = 0.06 m/s shows the highest volume fraction, whereas the maximum at n = 350 rpm is shifted toward $V_{mag}$= 0.08 m/s. The volume fractions of the simulated and experimental velocity magnitudes align very well with the velocities of $V_{mag}$ = 0.45 m/s for n = 250 rpm and $V_{mag}$ = 0.60 m/s for n = 350 rpm, thereby covering the largest fraction of the bioreactor. The occurrence of higher velocity magnitudes steadily declines, reaching a maximum velocity magnitude of $V_{mag}$ = 0.47 m/s (simulation) and $V_{mag}$= 0.51 − 0.58 m/s (experimental) at $n$ = 250 rpm for less than 10⁻³% of the total volume. At $n$ = 350 rpm, the volume fraction depicting maximal velocity magnitudes of up to $V_{mag}$ = 0.66 m/s for the simulation and $V_{mag}$ = 0.70 − 0.80 m/s for the experiment is less than 10⁻³%. Deviations of the simulation and experiment in higher velocities may be due to the finer time-step resolution of the simulation, thereby averaging higher occurring velocities within the impeller discharge stream compared to the tip speed. Consequently, using the Eulerian approach to generate time-averaged velocity profiles is unsuitable to detect velocities larger than the tip speed. However, as these elevated velocities are present within the impeller discharge stream, only transient profiles or Lagrangian particle trajectories are capable of resolving the apparent velocities, which is the advantage of the approach later discussed in this paper.

Figure 6 shows the vertical velocity profiles of the experimental and numerical data. For four different radial distances (r = 20, 30, 40 and 50 mm) to the stirrer shaft, the time-averaged velocity components are averaged azimuthally and normalized to the respective stirrer tip velocity. The modified velocity profiles for the baffled and unbaffled systems show very good agreement for both stirrer frequencies.

However, in the radial discharged flow of the two impellers in the baffled configuration, the experimentally modified velocity components are somewhat underestimated. This effect may be due to an insufficient local and/or temporal convergence of the mean experimental velocity. Possible options for improvement may be either a longer sampling interval or a higher particle density, neither of which can currently be implemented due to the limited memory and the local resolution of the high-speed cameras. Nevertheless, the measured and simulated velocity profiles show very good agreement outside of the direct influence area of the stirrers. Since precisely, these areas account for the significantly larger proportion of the velocity components occurring in the system, it can be assumed that the simulated data are valid. Due to very strong optical depletion caused by the bottom geometry of the stirrer vessel, which could not be fully compensated by the geometric calibration, a direct comparison of the simulated and experimental data in the bottom region of the vessel is not feasible.

Validation of numerical simulations by time-averaged flow fields is a method commonly used, showing fair agreement with experimental data also in other studies [11,12,19]. To furthermore account for maximal obtained velocities and transient structures within stirred vessels, a validation of the numerical simulation on the basis of the Lagrangian particle trajectories will be shown in the following.

Flow structures behind the blades are visualized by the Q-criterion and displayed in Figure 7A. The simulation is capable of capturing the two, characteristic trailing-roll vortices that develop behind the upper and lower edge of each blade as described by Nienow and Wisdom [44] (see Figure 7B). Circular rotation behind the blades is indicated by the white streamlines. Considering a tip speed of $V_{tip}$ = 0.47 m/s, the highest velocities $V_{mag}$ > 0.68 m/s are not obtained directly at the impeller tip but within the lower circular vortex within the impeller discharge stream. Apparently, the lower circular vortex rotates in the direction of disk rotation, contributing to increased velocities, whereas rotation for the upper vortex is vice versa.

Figure 8A shows selected particle trajectories of the measured and simulated particle movements in the system. Only those particle trajectories are shown for which the following conditions are fulfilled. On the one hand, the trajectories must have been observed over a period of at least 1.25 s. On the other hand, the particles must have a certain location. This location is described by the area, which is spanned from the stirrer shaft to the inner wall of the reactor at the height of the upper Rushton. The height of the surface is $h$= 12 mm. In addition, the area is oriented orthogonally to the focal plane of the four cameras in the experiment since most particle trajectories were observed for this area. In total, 191 trajectories were identified in the experiment and 219 in the simulation, which fulfill the described conditions and are therefore eligible for further evaluation.

In Figure 8B, the respective velocity components (resolved in time) for the filtered trajectories of each particle are shown normalized to all occurring velocity as velocity density distribution. A class width dV_mag = 2.734 mm/s is chosen for the histogram representation, corresponding to a subdivision into 256 velocity classes. It is remarkable that the maximum of the most frequently occurring velocity components of about $V_{mag}$ = 0.08 m/s can be observed, both in the experiment and in the simulation. The chosen representation method has the decisive advantage that, in particular, the transient velocity components are included in the comparison and are not filtered as in the representation of the time-averaged velocity vectors. This becomes clear by the fact that in both the experiment and simulation, particle velocities occur which are clearly above the stirrer tip velocity of $V_{mag}$ = 0.47 m/s.

In addition, Figure 8B shows the absolute displacement of the selected particles within the observed period for the experimental and simulated particles as a fraction on all particle displacements. For this, the distance between the start and end positions of the respective particles is determined. A subdivision into 25 classes is chosen. Here, a small displacement can be attributed either to overall low mean velocities of the respective particles or to the dwelling of a particle within the coherent vortex structures. However, since the particle trajectories shown here had to have passed the area around the stirrer at least once in order to be taken into account for the evaluation, a small absolute displacement is more likely to be due to particles staying within a temporally stable vortex, which results in the particle approaching its initial position again within the observed time period. In contrast, the high values of displacement are due to particles being attracted either from above/below the Rushton stirrer and being transported below/above. Again, the agreement between the experimental and simulated data is remarkably good.

In contrast, Figure 9A shows the particle trajectories, which were filtered analogously to those in Figure 8 and evaluated with respect to the velocity components and displacements but observed for the case without baffles.

For the setup without baffles, 308 as well as 376 particles fulfilling the described filter conditions can be identified in the experiment and in the simulation. The increase in the number of observed particles compared to the case with baffles is probably due to the stronger azimuthal flow, which significantly increases the probability that a particle migrates through the filter plane. This becomes evident when analyzing the trajectory lines in Figure 9A, where all path lines follow the azimuthal flow on a circular path around the stirrer shaft. Likewise, the maximum occurring velocity components of the particles are lower than in the setup with baffles, but still higher than the stirrer tip velocity. Compared to the data in the system with baffles, however, a small offset between the most frequently occurring velocity components can be observed at this point. In the experiment, the maximum of the velocity density distribution is at about $V_{mag}$ = 0.11 m/s, whereas the maximum of the velocity density distribution in the simulated particle trajectories is at about $V_{mag}$ = 0.08 m/s. This effect could also be observed for the time-averaged velocity fields. Nevertheless, there is a good similarity between the two velocity density distributions.

Furthermore, there is very good agreement between the simulated and experimental particle motions with respect to the absolute displacement, wherein, at this point, small absolute displacements are not due to dwelling within small temporal vortices. Instead, those particles have rotated once around the stirrer shaft on a certain radius without moving over the height or the radius, whereas longer absolute displacements are due to an exchange between the upper and lower compartments around the stirrer.

4. Conclusions

A 3 L laboratory vessel was examined by numerical simulations, and the numerical results were validated by experimental data. Power number determination by steady simulations with the RANS k − ε model MRF showed fair agreement with less than 11% deviation, on average, compared to experimental torque measurements. If computational power is limited, this approach reflects a computationally inexpensive method for equipment characterization, being fully automatable within the simulation software ANSYS Fluent. Nevertheless, in order to detect transient flow structures, necessary for reliable mixing time and gas transfer determination as well as to detect maximal occurring velocities and stresses, a time-resolved numerical scheme is mandatory, which was realized by an LB LES approach within the simulation software M-Star CFD in this study. Power number determination showed higher accuracy when averaging over different grid sizes; however, deviations for the chosen grid with reasonable computational demand accounted for less than 10% compared to the experimental data. As LB LES are not yet widely used for stirred vessel characterization in the pharmaceutical industry, time-averaged flow fields of the laboratory vessel were compared to 4D PTV data, showing very good agreement in the overall flow predictions for different reactor setups (baffled and unbaffled) and agitations frequencies. Additionally, LB LES in combination with 4D PTV data offer the unique possibility to validate transient structures by analysis of Lagrangian particle trajectories. To our knowledge, such a validation approach has not been published before. Both the distribution of the velocity magnitude and absolute displacement of the tracked particles show good agreement between the LB LES simulation and the experimental 4D PTV data. Furthermore, LB LES simulations are shown to be applicable for reliable and detailed prediction of mixing times [28]. Consequently, the validity of the obtained results justifies the transfer of the method to large, industrial scales.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/pr9060950/s1, Figure S1: Polyhedral mesh generated by Ansys Fluent meshing. From left to right, increasing mesh density of 100,000, 400,000 to 700,000 polyhedral cells Figure S2: Time- and azimuthal-averaged profiles of normalized velocity magnitude for the baffled system at $n$ = 250 rpm, Figure S3: Time- and azimuthal-averaged profiles of normalized velocity magnitude for the baffled system at n = 350 rpm, Figure S4: Time- and azimuthal-averaged profiles of normalized velocity magnitude for the baffled system at $n$ = 350 rpm, Figure S5: Time-averaged velocity magnitude vector field of the baffled system at $n$ = 350 rpm (A) and the unbaffled system at $n$ = 250 rpm (B).

Author Contributions

Conceptualization, M.K., J.F., T.W. and M.S.; methodology, M.K., A.v.K. and J.F.; software, M.K.; validation, J.F. and M.K.; formal analysis, M.K., A.v.K. and J.F.; investigation, M.K., A.v.K. and J.F.; resources, T.W. and M.S.; data curation, M.K. and J.F.; writing—original draft preparation, M.K. and J.F.; writing—review and editing, A.v.K., M.H., M.S. and T.W.; visualization, M.K. and J.F.; supervision, M.S. and T.W.; project administration, M.S. and T.W. Parts of this manuscript have already been presented as a plenary talk by J.F. and M.K. at the 2021 annual meeting of the ProcessNet Computational Fluid Dynamics and Multiphase Flow groups. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

The authors gratefully acknowledge Johannes Wutz for his support with numerical simulations using M-Star.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviation

Latin		Greek
$C$	Off-bottom clearance, m	$\Delta C$	Impeller spacing, m
$Co$	Courant number	$\Delta t$	Time step, s
$d$	Impeller diameter, m	$\Delta x$	Grid spacing, m
$D$	Tank diameter, m	$\eta$	Dynamic viscosity, Pa s
$f$	Probability density function	$\rho$	Density, kg m−3
$g$	Gravitational acceleration, m s−2	$\mathsf{\Omega }$	Collision operator
$H$	Tank height, m	$\epsilon$	Energy dissipation rate, m2 s−3
$h$	Surface height, m	$\sigma$	Surface tension, N m−1
$M$	Torque, N m
$n$	Agitation rate, rpm
$p$	Pressure, Pa
$P$	Power, W m−3
$P$ 0	Power number
$P_t$	Power by torque, W
$P_{\epsilon }$	Power by energy dissipation, W
r	Radial distance, m
Re	Reynolds number
$s$	Spatial resolution, m
$t$	Time, s
$U$	Three-dimensional velocity vector, m s−1
$V_{tip}$	Tip speed, m s−1
$V_{mag}$	Velocity magnitude, m s−1
$V$	Volume, m3

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Geometry of the bioreactor. Shaft mount (a), lid (b), bearing box (c), baffles (d), shaft (e), impeller (f), vessel (g).

Figure 2. Mesh study and variation of rotating reference frame volume. Evolution of power number by torque over number of grid cells for a fixed rotating reference frame volume (A) and depending on the ratio of impeller diameter or height to diameter or height of the RRF volume respectively (B) with the experimental value of Po = 7.5, shown as dashed line.

Figure 3. Time and azimuthal averaged profiles of normalized velocity magnitude. Profiles are averaged at different radial positions within the impeller discharge stream of top and bottom impeller.

Figure 4. Evolution of power number Po determined by torque Pt and by integral predicted energy dissipation Pε over number of grid points (A) with the experimental value of Po = 7.5 as dashed line (simulated average Po= 7.5 ± 0.5) and Po as function of Re for the baffled system (B).

Figure 5. Comparison of simulated (left of the stirrer shaft, mirrored) and experimental (right of the stirrer shaft) velocity field depicted as slice between two baffles of the baffled system at n = 250 rpm (A). Experimental (red) and simulated (blue) velocity distribution (3D) of the baffled system for n = 250 rpm (crosses) n = 350 rpm (circles) (B).

Figure 6. Comparison of the dimensionless axial time- and azimuthal-averaged velocity profiles.

Figure 7. Time-averaged (at fixed blade to baffle position) trailing vortices simulated with M-Star (A) and adopted from Van’t Riet [45] (B). Reprinted with permission from ref. [45]. Copyright 2021 Elsevier.

Figure 8. Visualization of experimental (magenta) and simulated particle trajectories for the baffled system at n = 250 rpm (A). Velocity distribution (upper) and absolute displacement of filtered particle trajectories (B).

Figure 9. Visualization of experimental (magenta) and simulated particle trajectories for the unbaffled system at 250 rpm (A). Velocity distribution (upper) and absolute displacement of filtered particle trajectories (B).