Industrial Case-Study-Based Computational Fluid Dynamic (CFD) Modeling of Stirred and Aerated Bioreactors.

Industrial bioreactors featuring inadequate geometry and operating conditions may depress the effectiveness and the efficiency of the hosted bioprocess. Computational fluid dynamics (CFD) can be used to find a suitable operating match between the target bioprocess and the available bioreactor. The aim of this work is to investigate the feasibility of addressing bioreactor improvement problems in the bioprocess industry with the aid of such mainstream tools as industry-standard CFD. This study illustrates how to effectively simulate both the impeller rotation and air supply and discusses the way toward model validation at the 4.1 m3 capacity scale. Referring to experimentally measured process values, the developed full-scale model successfully predicted the power draw, liquid phase level, and mixing time with errors lower than 4.6, 1.1, and 6.7%, respectively, thus suggesting the illustrated approach as a best practice design method for the bioprocess industry. The validated model was employed to improve performance by reducing the power draw in aerated conditions with a minimal operational derating.

Introduction

Bioreactors denote equipment designed to provide biological, biochemical, and biomechanical conditions that satisfy the needs for the optimal growth of biomass (i.e., animal or plant cells, tissues) in terms of quantity and quality and/or for carrying out biochemical reactions on an industrial scale suitable for the synthesis of desired products.[1] Regarding microalgal cultivation, fermentation represents a long-established type of a highly controlled production process aimed at optimizing cell growth and the production of high-value-added products.[2−4] The closed structure of fermenters makes it easier to maintain axenic conditions that are mandatory to ensure the product safety that is required by the food and other industries; furthermore, the high cell densities that can be achieved in fermenters (>100 g L–1) put this technology in the front line for industrial-scale microalgal culturing processes.[5] Mechanically stirred aerated bioreactors are widely used in the production of metabolites from cultures under aerobic conditions, in which the product formation rate is mainly dictated by the gas–liquid mass transfer rate. Indeed, the dissolved oxygen concentration is the main parameter that should be maintained at the value required by the specific process to gain the maximum performance from the cultured aerobic microorganism. Usually, air is introduced from the bottom by injection through a single nozzle or one or more spargers. The bubbles generated from the sparger are entrained in the circulating liquid flow, and by managing the stirring rate, they are recirculated through the stirrer, thus increasing both the gas holdup and the contact time between the liquid and gas phases and warranting the required volumetric overall oxygen transfer coefficient kLa.[6,7] The mechanical stirring system consists of a single impeller or a combination of several impellers of the same type (triggering axial or radial flows) or of different types whose arrangement and installation address the needs of the specific culture process and ensure the required mixing of the culture.[8]

Frequently, an installed mechanically agitated fermenter needs to be adapted to some fermentation process that needs to be scaled up.[9,10] Ideally, the existing equipment will be used at its best without any limitations, but chances are that some modifications are required either to operational variables (i.e., by increasing the oxygen concentration in the air supply) or to some mechanical parts (i.e., replacement of the agitator impeller to improve mixing in some aspects). Predicting the optimal performance that can be achieved in a bioreactor requires a thorough knowledge of its hydrodynamics;[11,12] otherwise, one can end up in a lengthy and costly trial-and-error work. This paper presents the use of computational fluid dynamics (CFD) to model the hydrodynamics of bioreactors at the industrial (>1 m3) scale to analyze and improve the extent and uniformity of the mass transfer coefficient. The adopted case study is an aerobic fermenter with mechanical agitation aimed at culturing microalgae. It was developed within the framework of the MEWLIFE European project and is installed in the NextChem premises in Rome.[13] Two rotation speeds were investigated during the model development and validation phase of the work, namely, 116 rpm, which is the design rotation speed, and 58 rpm. A further one (87 rpm) was adopted to demonstrate the role of CDF in improving the bioreactor management. A Eulerian–Eulerian multiphase model was used to describe the two-phase flow. The k-ω shear stress transport turbulence (SST) model was employed to characterize the fluid dynamic features of the reactor,[14] which is different from other studies that employ the standard k-ε model.[15,16] The k-ω SST turbulence model was chosen for closing the Reynolds-averaged Navier–Stokes equations (RANS) since it has shown a remarkable robustness and reliability in simulations involving complex flows.[17]

An extensive experimental campaign was carried out to measure key parameters such as the power draw, liquid phase level, and mixing time to compare to numerical outcomes and to obtain a validated numerical model. It is worth noting that the works available in the literature estimate the behavior of industrial-scale reactors only theoretically;[18] focus only on one critical aspect at a time, generally either mixing or aeration;[19,20] or deal with different reactor geometries,[21−24] mostly at the bench or pilot scale. So far, there are no experimentally validated models in the literature dealing with industrial-scale bioreactors encompassing both air dispersion and mixing of the aerated mixture. Due to the high computational time required for a single complete simulation, only small-scale reactors are typically modeled.[25−29] When model validation at the bench scale is carried out, researchers are free to choose whatever instruments fit their needs to ensure a thorough characterization of all the operating variables that can be useful for model validation, even beyond what is normally required for fermentation control (or can be installed and operated) in production units. Differently, model validation against an installed industrial site bioreactor may often be hampered by constraints determined by the installed equipment or instruments. The availability and location of free nozzles, the specific features and installation location of the installed instruments, and the large volumes may require that some aspect of the validation procedure be adapted to the adopted industrial unit. Furthermore, such accurate analysis tools as the particle image velocimetry[30] cannot be applied due to the lack of optical access through the metal tank shell.

The validated CFD model was finally employed to improve the performance of the fermenter by searching for a compromise between satisfactory uniformly distributed values of kLa, identified as a critical parameter for the process, and acceptable costs of management of the bioreactor itself. This was achieved by investigating both the two rotation speeds used for the validation procedure and an intermediate rotation speed, i.e., 87 rpm.

The main innovation brought by this work lies in demonstrating the usability of CFD as a suitable tool for dealing with large-scale bioreactors and the improvement of their operation. This purpose is achieved by fully describing the extensive validation procedure carried out and illustrating the pathway to obtain industrial-scale relevance with out-of-the box CFD tools.

Material and Methods

Case Study

The present work is based on an ″equipment rating″ approach and adopts an industrial-scale mechanically agitated aerobic fermenter as the case study. The equipment is devoted to the heterotrophic cultivation of microalgae and is part of a pilot plant designed and used in the framework of the (currently ongoing) MEWLIFE European research project (Figure ).

Figure 1

Bioreactor located at the MEWLIFE plant in the NextChem site in Rome.

It features a cylindrical shell with two elliptic bottoms and a total (geometric) volume of 4.1 m3. Internally, the vessel presents three longitudinal baffles with a 120° radial arrangement installed next to the cylindrical shell, a toroidal gas sparging system with downward holes, and a steel pipe for the level measurement system. The agitation system includes an axial shaft on which an upper Rushton turbine and a lower downward-pumping pitched-blade turbine are installed. The design-time motivation for the bottom location of the axial impeller was to prevent solids settling on the vessel bottom while ensuring an even liquid flow throughout the lower and bottom sections of the vessel. The choice of type and size of the impellers, and their location had been made at design and were not available as adjustable items during the time frame of the present research. The bioreactor is schematically shown in Figure , where its relevant parts and their size are also reported.

Figure 2

Overall configuration of the bioreactor implemented as a CAD model: (a) perspective, (b) front, and (c) from above. Details of the stirring system: (d) Rushton turbine and (e) impeller with four inclined blades.

Numerical Model Implementation

The core of a simulation carried out via computational fluid dynamics lies in the solution of the Navier–Stokes equations. One of the possible solutions of the Navier–Stokes equations is the Reynolds-averaged Navier–Stokes (RANS) equations, where the instantaneous transport properties are considered as the sum of their mean value (time-averaged) and a fluctuating quantity due to turbulence. Furthermore, the presence of a nonlinear Reynolds stress term requires additional modeling to close the RANS equations, and this has led to the creation of many different turbulence models. According to the Boussinesq hypothesis (turbulent-viscosity hypothesis),[31] the Reynolds stresses can be correlated to the mean velocity:where ρ is the density, k is the turbulence kinetic energy, δ is the Kronecker delta, μt is the turbulence eddy viscosity, and v̲i the mean velocity component along direction i.

In this work, the simulations were conducted in transient mode, choosing the following models:

Volume of fluid (VOF) with the implicit body force option for modeling the interaction between the phases (air and water in this case), VOF being a reliable model with a reduced computational cost. This Eulerian model considers the phases as continuous and non-interpenetrating. The model performs the calculation of one continuous phase, calculating the other one as a difference. This calculation approach guarantees a low computational cost under the restrictive hypothesis, verified in the present case unit, that the dispersed phase is less than 10% of the total.[32,33] Gas holdup in stirred fermenters can range from 1 to 20%.[34] Therefore, the proposed approach is suitable for average gas rate sparged units.

k-ω Menter’s shear stress transport (SST) for turbulence modeling, as it combines a robust formulation and near-wall region accuracy (imported by the k-ω model) with an independence from free flow in the bulk areas (guaranteed by the k-ε model). Furthermore, it is also suitable for a wide range of Reynolds numbers (even low ones). Finally, it also guarantees a low computational cost.[35,17]

Multireference frame (MRF) for modeling the rotation of the impellers. MRF prescribes that the equations are expressed in a reference frame that rotates at the same rotation speed of the impeller and are solved for a stationary mesh; this implies that the technique has a computational cost lower than the sliding mesh approach, where rotation is simulated by means of a rotating grid.[35−38]

The ANSYS Fluent 2019 R1 software was used to perform the simulations, and the bioreactor geometry was reconstructed in CAD using the SpaceClaim software within ANSYS (Figure ).

Boundary Conditions

Boundary conditions are required along all the boundaries of the solution domain to define the specific fluid flow. A ″Wall″ condition was set on the walls of the bioreactor and on the solid surfaces of all components, such as the impellers, shaft, sparger, baffles, bottoms, still pipe, pipes that feed acids and alkali, and bolted joint. The ″Velocity Inlet″ condition was set on the holes of the sparger, setting the air fraction volume equal to 1. Atmospheric pressure, via the ″Pressure Outlet″ condition, was specified at the dead end of the still pipe and at the bioreactor top. Along these latter two surfaces, a backflow volume fraction of air only was imposed. Finally, the boundary condition named ″internal″ (which serves to delimit the zone where the flow rotation takes place) was set at the walls of the rotation zones modeled with MRF.

Mesh Properties and Grid Validation

The volumetric mesh was created by means of Fluent Meshing using the Poly-Hexcore function. Solutions over different grid resolutions were computed and compared to determine the optimal cell number and size. The final grid consisted of approximately 940,000 cells, a number that was considered an acceptable compromise between calculation accuracy and computation times. In fact, a further grid refinement was considered by decreasing the size of cells by 20%. The number of cells of the refined grid was approximately 1,100,000. For the 116 rpm case, the power draw, liquid phase level, and mixing time were computed with both grids and compared; differences happened to be lower than 2.3, 1.0, and 3.2%, respectively, in front of an increase of the computation time of roughly 15%.

The main characteristics of the mesh that play a major role in the accuracy and stability of the numerical computation are reported in Table .

Table 1

Main Characteristics of the Created Mesh

minimum volume [mm3]	0.26
maximum volume [mm3]	3.65 × 104
maximum skewness	0.8
minimum orthogonal quality	0.4

For 3D simulations, a maximum skewness less than 0.9, as well as a minimum orthogonal quality value around 0.4, certifies a good quality mesh.

Rotating Zone Selection

To implement the multiple reference frame approach, a region of cylindrical shape containing each impeller has to be identified. For this purpose, the internal volume of the bioreactor was divided into three zones: two rotation zones (one for each impeller) and one for the remaining volume of the reactor. The sizes of these rotation zones were carefully chosen. In fact, by increasing their size, the mixing effect of the fluid inside the reactor increases. The ratio of the diameter to the height of the impeller and that of the diameter to the height of the rotation zone were defined based on previous studies in the literature and, in particular, on the work of Patil et al.[37] Patil et al. investigated the effect of adopting rotation zones of different sizes for an axial and a radial impeller in a stirred vessel of the same shape and volume as the reactor under examination. As a result of this study, a rotation area with diameter da = 1.182 m and height ha = 0.186 m for the axial impeller and a rotation area with diameter dr = 0.945 m and hr = 0.246 m for the radial impeller were created (Figure ). The comparison between quantities measured via the experiments and calculated by the numerical model represents an indirect approach to verify the consistency of such a choice.

Figure 3

Rotation zones’ position and sizes.

Experiments for Model Validation

Specific experiments were designed and carried out to validate the described CFD model. The quantities adopted for the validation are the (i) stirring power requirement, (ii) mixing time, and (iii) level of the liquid phase inside the bioreactor in the presence of air bubbles, i.e., gas holdup.

Stirring Power Requirement

The calculation of the stirring power requirement with the numerical model is straightforward and easily comparable with the values obtained through experimental measures. As Ascanio et al.[39] and Gogate et al.[40] report, the most frequently used experimental techniques for the evaluation of energy consumption in agitated tanks and fermenters employ watt meters, ammeters, colorimeters, dynamometers, torsiometers, and systems based on strain gauges. Each of these systems has its own advantages and disadvantages and may be chosen based on factors such as investment costs, accuracy, and measurement range. In this work, the effective method proposed by Taghavi et al.,[41] which compares the power values estimated by CFD simulation with the electric power absorbed by a gear motor provided by the inverter controller, was adopted. Electrical and mechanical losses were considered. Specifically, the experimental measurements for determining the power required for mixing consisted of the measurement of the intensity of the current absorbed by the gear motor (BOX110 from Motive) that was supplied by the inverter (INVT GD-20 EU) installed for the rotation speed control. Using the voltage applied V and the current phase shift cos φ declared by the engine manufacturer (0.85), the absorbed power was calculated using the following equation from basic electrical engineering textbooks:where I stands for the absorbed current. To calculate the power required for stirring the fluid only, experimental tests with the impellers rotating in the empty reactor were carried out at the same rotation rate. The power measured with those tests was assumed to be equal to the sum of electric and mechanical losses, so these values were subtracted from those obtained in the experimental measurements with the full reactor. In this way, the experimental values can be compared to the numerical predictions.

The following experimental tests were repeated three times to verify the repeatability of the outcomes:

The time history of the current absorbed by the electric engine in the empty bioreactor was recorded from the onset of the experiment until the stationary conditions were reached for both the 116 and 58 rpm rotation speeds;

The time history of the current absorbed by the electric engine in the bioreactor filled with water (3 m3) from the onset of the experiment to the rotation steady-state condition was recorded for both the 116 and 58 rpm rotational speeds; additional experiments were carried out with air supply.

The model prediction of power draw was calculated (in watts) as the product of the moment calculated (by Fluent) times the corresponding rotation rate by applying any needed conversion.

Mixing Time

Conductivity tests were carried out to determine the mixing time within the bioreactor for the two different stirring rates. A concentrated saline solution was rapidly injected in the reactor just below the free surface, and the electrical conductivity (proportional to the raw resistance measured) was recorded over time after the injection.[19,20,29] The experimental tests were carried out by using a conductivity meter (Mettler Toledo Cond 7100e) that permitted reading the instantaneous resistance without applying any digital filter and with a negligible measurement delay (<1 s based on the performed checks). The validation tests were carried out twice for each rotation speed, each with a different volume of a concentrated saline solution (300 g/L). The time history of the measured electrical resistance was recorded after the injection of the saline aqueous solution through a nozzle at the top of the fermenter. The visual inspection of the time history of conductivity and the analysis of the relevant data were used to determine the values of circulation and mixing times.

The numerical evaluation of mixing time was made by a Lagrangian approach (circulation and mixing times were calculated using the particle trace simulation technique, whereby the trajectory of a set of nonreactive massless virtual particles is analyzed). After the steady condition of the system relevant to stirring was reached (40 s was assumed as a safety value for both the rotation speeds), the time at which the trajectories of the particles evenly occupy the entire volume of the bioreactor was identified. The particle motion was simulated in ″transient mode″ selecting the ″Interaction With The Continuous Phase″ option, and both the 116 and 58 rpm rotation speeds were investigated. The application of this approach was used by Villiger et al.[42] to trace the pH profile in different volumes of a fluid in motion; Villiger et al. concluded that such use of the Lagrangian method leads to results almost identical to the Eulerian one used in earlier works.[19,20]

Level of the Liquid Phase

To validate the gas holdup prediction, several experimental measurements were conducted using a radar level gauge (Micropilot FMR51 HART, Endress+Hauser, Switzerland). First, the bioreactor was brought to a steady state in terms of stirring. Then, air insufflation was started, and the water level was recorded until a new steady state was attained. The measurements were carried out for both rotation speeds (116 and 58 rpm).

From the numerical point of view, the interface between water and air at steady state of each simulation does not appear sharp due to the intrinsic limit of the Eulerian–Eulerian model like the volume of fluid employed to simulate the two-phase flow. As Wang et al.[43] found out, to make the numerical results comparable to the experimental ones when using a Eulerian–Eulerian model, it is mandatory to set a specific value for the air volume fraction to identify the interface position. Before addressing the gassed condition, the average liquid phase level in steady-state conditions was evaluated as a spatial average using six planes passing through the central axis of the fermenter. Then, following Wang et al., the interface was detected ranging the air volume fraction values between 0.4 and 0.6 and comparing the calculated position of the interface to that provided by the experiments.

kLa

The capability of providing an adequate oxygen supply is a fundamental aspect in designing a bioreactor aimed at aerobic fermentations. The mass transfer coefficient in the liquid phase, kL, is one of the fundamental parameters to be analyzed to verify if the entire liquid volume is aerated to a sufficient extent for the design biomass concentration and with a sufficient uniformity so that local limitations will not be reached when the concentration approaches the design value.

It is known that the kLa value may vary within the bioreactor, and consequently, it is important to consider not only its average value but also its local value. The kLa value was evaluated as an outcome of the numerical model by separately calculating the local mass transfer coefficient, kL, and the local specific surface area, a, and taking their product as the local mass transfer coefficient. The local mass transfer coefficient was calculated by using the theory of penetration:where D stands for the diffusion coefficient of oxygen in water (which at room temperature is equal to 100 × 10–11 m2/s), ρL and μL are the density and viscosity of water (1000 kg/m3 and 8.9 × 10–4 Pa s3, respectively), and ε is the turbulent kinetic energy dissipation rate computed by CFD. The local specific contact surface between the gas and the liquid was computed aswhere φ is the local volume fraction of the gas and d is the average diameter of the bubbles, which was established to be equal to 4.4 mm. This value was calculated after ref (44) based on the orifice Reynolds number (11,700) and is consistent with ref (45) indicating a constant average diameter to be included in the equation, usually between 2 and 4 mm, when VOF is used to model the two-phase flow.

Results

Power Required for Agitation

The comparison between the values of absorbed power computed by the simulations and those measured with the experimental tests for the two case studies (116 and 58 rpm) both without and with air is reported below (Table ). A remarkable agreement is observed between the numerical and experimental data at steady-state conditions for all the four cases investigated.

Table 2

Mean Absolute Percentage Error for Absorbed Power

case	116 rpm no aeration	58 rpm no aeration	116 rpm aeration on	58 rpm aeration on
mean and standard deviation (experimental values)	1076 ± 25 W	153 ± 4 W	855 ± 27 W	107 ± 4 W
mean and standard deviation (numerical values)	1071 ± 38 W	153 ± 6 W	873 ± 5 W	112 ± 6 W
mean absolute percentage error (MAPE)	1.0%	0.9%	2.5%	4.6%

To evaluate the quality of this agreement, the mean absolute percentage error (MAPE) was computed employing the experimental and the computed values at steady state (after 40 s for the case of stirring only and after 60 s for the case of stirring with air insufflation). The MAPE values presented in Table confirms the high accuracy of the CFD model to reproduce the experimental results for power consumption.

Mixing Time Results

Figure displays the experimental results obtained by measuring with the conductivity probe the variation of the electrical resistance after the addition of 500 mL of an aqueous NaCl solution at 300 g/L. For both tests, the profiles obtained show oscillating trends, which allow us to obtain information about the mixing time useful for our analysis.

Figure 4

Measurement of electrical resistance over time following the addition of 500 mL of aqueous solution of NaCl at 300 g/L with a rotation speed of (a) 116 and (b) 58 rpm.

The two quantities used for validating the mixing features are circulation time and mixing time. When the time history of the measured resistance is considered, the circulation time is calculated as the time elapsed between any two successive peaks, while the mixing time is estimated from the time required for the overshoot beyond the asymptotic value of resistance to decay below 10% of the difference between the initial value of resistance and the steady-state one. In Figure a, referring to the rotation speed equal to 116 rpm, a circulation time and a mixing time of about 4.5 and 15 s were estimated, respectively. In Figure b, referring to 58 rpm, a characteristic time of about 9.5 s and a mixing time of about 24 s can be detected. It is important to underline that the weakness of the first peak on the plot of Figure b leads to the hypothesis that, at 58 rpm rotation speed, the full recirculation of the fluid cannot be reached. Evidence that confirm this hypothesis will be shown in a later section of this work.

To analyze the mixing effectiveness of the bioreaction environment and evaluate the mixing time with the numerical model, the simulation was carried out until the steady-state condition of the system for stirring only was reached, i.e., after 40 s from the start of the simulation. After that, 5000 particles were inserted within the flow domain using the discrete phase model (DPM) in transient mode, and the simulation was carried out for a further 40 s. The particles were introduced in a volume corresponding to the one employed for injecting the tracer during the experimental tests.

A preliminary qualitative analysis of the trajectories was carried out that consisted of displaying them at different residence times and identifying the time at which the trajectories of the particles fully occupy the entire volume of the bioreactor. For both rotation speeds, the particles compactly follow the flow dictated by the agitation system, and after a certain time, they start to follow trajectories that differ from one another and eventually distribute evenly over the vertical direction, i.e., over the entire liquid volume. Starting from the assumption that mixing can be considered ″complete″ when particles trajectories ″fill″ the entire bioreactor volume, a mixing time of roughly 15 s for the design rotation speed (i.e., 116 rpm) and 25–30 s for half the design rotation case was detected. These results are remarkably in accordance with the results of the corresponding experimental tests (Figure ). To further corroborate these results, a quantitative analysis was conducted to determine the average distribution of the particles along the vertical direction as a function of time. The mixing times derived from this analysis are approximately 16 s for the case at 116 rpm and about 23 s for the case at 58 rpm, values that are in accordance with the experimental values. Furthermore, by observing the initial part of the diagram, the concept of circulation time was generalized to the average time required for particles to return to the same vertical position. Thus, the circulation time in the case of the 116 rpm rotation speed was calculated to be equal to 4.6 s, and in the case of the half rotation speed, it was calculated to be 7.2 s (Figure ). Again, these values are coherent with those estimated by the experimental validation test reported in Figure . To quantify the correspondence between experimental and numerical values, Table summarizes the results obtained and presents the outcome of the MAPE analysis. Even in this case, the remarkable correspondence between model results to experimental data is shown. The only comparison that leads to a high mean absolute percentage error is the circulation time comparison at 58 rpm, which appears to validate the hypothesis that, in this condition of rotational speed, a perfect recirculation was not reached.

Figure 5

The z coordinate of the particles as a function of time for (a) 116 and (b) 58 rpm.

Table 3

Mean Absolute Percentage Error for Mixing Time and Circulation Time

case	116 rpm	58 rpm
mixing time (experimental values)	15 s	24 s
circulation time (experimental values)	4.5 s	9.5 s
mixing time (numerical values)	16 s	23 s
circulation time (numerical values)	4.6 s	7.2 s
mixing time mean absolute percentage error (MAPE)	6.7%	4.2%
circulation time mean absolute percentage error (MAPE)	2.2%	24.2%

Liquid Phase Level

Table shows the results of the liquid phase level analysis once the steady state was reached. This analysis made it possible to confirm that 0.5 is a reasonable value for the water fraction (wf) value to identify the gas–liquid interface considering that the MAPE values of the liquid phase lever value obtained with wf = 0.5 are the lowest ones.

Table 4

Mean Absolute Percentage Error for the Liquid Phase Level

case	116 rpm	58 rpm
liquid phase level (experimental values)	1.763 m	1.708 m
liquid phase level wf = 0.4 (numerical values)	1.822 m	1.742 m
liquid phase level wf = 0.5 (numerical values)	1.782 m	1.690 m
liquid phase level wf = 0.6 (numerical values)	1.700 m	1.637 m
mean absolute percentage error (MAPE) with wf = 0.5	1.1%	1.0%

Bioreactor Performance Improvement: kLa

After model validation, the bioreactor performance was investigated to identify areas of potential improvements. The volumetric overall oxygen transfer coefficient kLa was estimated by considering its spatial distribution after an averaging procedure along six planes passing through the axis of the bioreactor was performed and considering five different times, all pertaining to the steady-state condition. A time average was performed because a single instant would have provided a snapshot of the position of the bubbles within the bioreactor volume. Averaging over several times allowed the characterization of all parts of the control volume and made it possible to verify that they were well oxygenated.

Figure shows the local distribution of kLa values obtained for the investigated rotation speeds, i.e., 58 (Figure a), 87 (Figure b), and 116 rpm (Figure c). Table presents three bulk average values computed considering three characteristic regions, i.e., the area above the Rushton turbine (area #1), the area between the two impellers (area #2), and the area below the impeller with four inclined blades (area #3). According to refs (46−48), kLa values around 0.02–0.03 s–1 are optimal for the growth and maintenance of a microalgae culture.

Figure 6

kLa maps (s–1) for the rotation speeds (a) 58, (b) 87, and (c) 116 rpm.

Table 5

Bulk kLa Values in Three Characteristic Regions within the Fermenter

	116 rpm	87 rpm	58 rpm
area above the Rushton turbine (area #1, 52% of the total volume)	0.018 s–1	0.021 s–1	0.019 s–1
area between the two impellers (area #2, 35% of the total volume)	0.019 s–1	0.015 s–1	0.010 s–1
area below the impeller with four inclined blades (area #3, 13% of the total volume)	0.014 s–1	0.010 s–1	0.003 s–1
overall kLa value	0.018 s–1	0.018 s–1	0.014 s–1

The correct operation of stirred aerobic fermenters requires that the stirring rate is well beyond the value for which the impeller is flooded at the prevailing gas rate and that mass transfer covers the expected demand of the microbial culture. The transition between ″flooding″ and ″loading″ and from ″loading″ to ″complete gas dispersion″ is described in terms of the impeller gas number Fg (Fg = Qg/nD3) and of the Froude number Fr (Fr = n2D/g) by the following general relationship:where A, B, and C assume the set of values (30, 3.5, 1) for the transition from flooding to loading and (0.2, 0.5, 0.5) for the transition from loading to the complete dispersion of the gas. Based on these two equations, the minimum stirring speed that warrants that flooding is avoided at the prevailing gas rate (0.96 vvm) is 16 rpm, while 127 rpm is predicted to be required for a complete dispersion of the gas.[49] Therefore, the rated unit operates at nearly complete dispersion of the gas. Indeed, the CFD predicted that the kLa at 50% of the design speed (the flooding–loading transition value) exhibits a significantly non-uniform distribution over the bioreactor volume, with ample zones where the local estimated value of kLa is markedly lower than the average calculated value (dark areas in Figure a plot). Indeed, the coefficient of variation of kLa is 43%. Both experimental and predicted mixing and circulation profiles demonstrated the inability of the impellers to warrant a thorough circulation. Contrary to this, the design rate kLa exhibits a fairly good homogeneity (Figure c), with a coefficient of variation less than 10% and a sensibly higher average value (by about 30% of the value predicted at the lower stirring rate).

It should be noted that mass transfer limitation, even if it is confined to a specific zone of the bioreactor, is destined to eventually create a bottleneck for the production deployment of the rated equipment when the biomass load is increased to maximize the per-batch production capacity. A bottleneck is eventually met because the maximum volumetric concentration of microbial biomass that can be sustained in fully aerobic metabolism meets an upper ″cap″ dictated by oxygen starvation. Indeed, the culture medium is fully saturated with oxygen only in the vicinity of the impellers (where the shear is most intense), and even there, the oxygen concentration that can be attained is modest (as an effect of the scarce solubility of oxygen in water and its further depression due to the culture medium temperature and composition). It should be considered that the microbial biomass is entrained by the liquid and circulates together with it, and by consuming oxygen at its characteristic rate (qO2, in kgO2 per kgbiomass and unit time), it lowers the local oxygen concentration of the suspending liquid volume element it is traveling in. Thus, after having been brought to the saturation value near the impellers, the oxygen concentration ″sensed″ by the microbial biomass decreases during its travel across the bioreactor volume until it returns to the high shear zones. By also considering that cell metabolism changes (this change may be anywhere from a slowdown to a more pronounced shift) below the so-called ″critical″ oxygen concentration value, it can be readily concluded that an upper ″cap″ of biomass productivity exists for a given bioreactor volume-distribution of kLa and that it can be investigated by appropriately modeling the dissolved oxygen concentration decrease during the circulation of the suspension. By considering the oxygen mass balance, one may write

Equation can be solved with the initial condition cL = csat at time t = 0 when the liquid element leaves the high-shear (and also high-kLa) zone and the final condition cL = ccrit after one circulation time. By substituting the lowest predicted value of the kLa and then solving the balance for unknown X, one obtains Xcrit, that is, the maximum value of biomass concentration that the bioreactor can sustain at the operating conditions that warrant the substituted value of the kLa.

If the biomass concentration were to be increased beyond Xcrit, part of the reactor volume would perform below expectations, and other undesired consequences (metabolic shifts) might be triggered by oxygen scarcity.

Figure shows the maximum concentration of microbial biomass that the bioreactor can sustain under aerobic conditions requesting 0.28 gO2 per gram of dry microbial biomass per hour.[50] It should be noted that there is no such thing as a ″typical″ value of qO2, as the oxygen demand rate of a microorganism depends largely on the substrate it is growing on (the more reduced the substrate, the higher the oxygen demand rate), the growth rate, and the fact that biomass, or rather an exogenous metabolite, is the product. Since low-cost fermentation media are frequently sourced from vegetal materials such as starch, they are fairly oxygenated, and the qO2 value that has been adopted for the present scenario analysis of the adopted case study is appropriate. It can be predicted that at the lower stirrer rate, about 36 kg of dry biomass can be obtained at the end of the fermentation (12 kg/m3 × 3 m3), while this value increases up to 66 kg at the higher rate. Under these conditions, from 24 to 43% of the inlet oxygen is absorbed by the fermentation process and replaced by an almost equivalent amount of CO2. A further harvest increase can be obtained by increasing the oxygen transfer. This latter, in turn, can be warranted by increasing the driving force, the stirring rate, the gas supply, or a combination of these parameters. While the driving force can be increased by increasing the backpressure (which was assumed to be 1.2 bar in the present case) or by modifying the composition of the inlet air (e.g., by supplementing the air feed with pure oxygen), it should be noted that most microbial processes are carried out as batch processes, where the biomass concentration is initially very low (inoculum) and increases progressively (initially in a very slow manner and then exponentially) over time so that the limitation may only appear at the later stages of the fermentation when the biomass concentration approaches the maximum value. Therefore, any bioreactor management policy that entails a higher operating cost can be reserved for the ″final″ part of the batch process when the biomass concentration approaches the affordable limit. A defective mass transfer does not prevent a fermenter to be used but, by limiting the maximum sustainable biomass concentration, would limit the overall bioprocess profitability.

Figure 7

Maximum biomass concentration that can be sustained by the bioreactor in completely aerobic conditions in the upper, intermediate, and lower zones of the bioreactor as a function of the stirring rate and oxygen transferred per unit energy expended as a function of the stirring rate.

Figure also shows the mass transfer energy efficiency of the rated design. As can be seen, the amount of oxygen transferred per unit energy expended for the stirring is 0.232, 0.242, and 0.198 kg/kWh at the rotational speed of 58, 87, and 116 rpm, respectively. This number is in the order of magnitude usually found for similar equipment.

Conclusions

The aim of this work was to introduce a case-based feasibility study of CFD modeling for bioreactors at the industrial scale of holdup volume. A 4.1 m3 mechanically agitated bioreactor aimed at heterotrophic microalgae fermentation was adopted as the case study. The computational fluid dynamic model was implemented as a multiphase model using the Poly-Hexcore 3D gridding scheme, the volume-of-fluid model of interaction between phases, the k-ω model for turbulence, and the multiple reference frame model for rotating impellers in the commercial software ANSYS Fluent. The principal innovation brought by this work is validation in a ″real-life″ situation on a large-scale bioreactor using a rating approach and model-assisted improvement of operational mass transfer and biomass-sustaining capability of the case-study production unit. Cost-effective modeling and workable validation instrumental and model-based techniques were conceived and implemented.

The remarkable agreement between model and experimental results, centered on the evaluation of the mixing effectiveness (circulation and mixing time) and of the power requirement (power draw of the stirrer) of the modeled unit, demonstrates the validity and the effectiveness of the adopted approach. Indeed, the power draw, liquid phase level, and mixing time had mean absolute percentage errors lower than 4.6, 1.1, and 6.7%, respectively.

The results of the simulations afforded a valid compromise between an optimal distribution of satisfactory values of kLa and sustainable operating costs: at 75% of the design rotation speed, a uniform distribution of the mass transfer coefficient near the optimal range was obtained, along with a significant operation energy saving (−50% with respect to the design speed and the optimal mass transfer per total expended energy performance), at a modest production derating (−14%).

The present work paves the way to the full industrial-scale bioreactor design and performance improvement by CFD modeling based on commercial CFD software.