Numerical Simulation of Bubble-Liquid Two-Phase Turbulent Flows in Shallow Bioreactor

: An improved second-order moment bubble-liquid two-phase turbulent model is developed to predict the hydrodynamic characteristics of the shallow bioreactor using two height-to-diameter ratios of H / D = 1.4 and H / D = 2.9. The two-phase hydrodynamic parameters, the bubble normal and shear stress, the bubble energy dissipation rate, the bubble turbulent kinetic energy, etc. were numerically simulated. These parameters increased along with ﬂow direction and constituted a threat to cells living at far distance away from the gas jetting inlet regions, rather than a ﬁnding of higher cell damage at near the jetting inlet region, as reported by Babosa et al. 2003. A new correlation named the turbulent energy production of bubble-liquid two-phase ﬂow was proposed to successfully verify this experimental observation. A smaller H / D ratio makes more contributions to the generation of lower turbulent energy productions, which are in favor of the alleviation of cell damage. The extremely long and narrow shape of the bioreactor is deteriorative for cell living.


Introduction
The shallow gas-liquid bioreactor characterized by the lower height-to-diameter ratios, has been widely applied in the field of biotechnology engineering, chemical engineering and pharmaceutical industries because of its low gas pressure drop, absence of moving parts, low cost and robust liquid phase residence time, especially eliminating the high-pressure drop of the tall reactors [1][2][3], Bubble-liquid hydrodynamics in the shallow reactor are predominantly controlled by sparger construction and the ratios of height to diameter (H/D), which are quietly different from those of a taller-shaped reactor. To thoroughly understand the hydrodynamic characteristics, it is essential to grasp the essences of optimization and scale-up strategies. However, the complex mechanism of bubble-liquid two-phase turbulent hydrodynamics and the effects of sparger design and H/D ratios on flows and transport characteristics have not been acquainted clearly. Up to now, most investigations have been performed heavily on an empirical or semi-empirical basis, and evidences that only originated from the experimental data, resulting in the difficulty to fulfill in-depth experimental validations under different platforms and derivations of theoretical applications due to their limitations [4][5][6].
With respect to the animal cell culture of biotechnology engineering, although early industrial exploitation began with Salk poliovirus vaccine production in primary monkey kidney cells, commercial application of animal cells to produce human biopharmaceuticals has rapidly increased in recent years [7,8]. Thus, reasonable design and optimized strategies of the bioreactor are the key issues for successfully producing commercial biopharmaceuticals. Both the experimental and the simulation investigations have been focused on the effect of the height-to-diameter (H/D) ratios on the flow performance in the shallow bubble column reactor. Throat et al. [9] investigated the combined effects As mentioned above, the standard k-ε model or the renormalization group RNG k-ε model are currently used in the ANSYS Fluent commercial software. Moreover, the existing two-fluid models are based on the isotropic turbulent viscosity. Meanwhile, liquid turbulence and bubble turbulence are distinctively anisotropic. Even though the EDR or the normal/shear stresses are considered as generally hydrodynamic parameters related to cell damage, the main findings are that they have larger values at a far distance away from gas jetting inlet regions. This conclusion failed to elaborate on the higher death rate. Furthermore, the bubble rupture or coalescence resulting in cell damage are also unreasonable because these events are primarily located at the top of the liquid surface or nearby the stirred blade, rather than at gas jetting inlet regions [29][30][31][32][33][34][35][36].
In this work, an improved second-order moment bubble-liquid two-phase model that can fully consider the anisotropic bubble dispersion characteristics was proposed to predict the bubble-liquid hydrodynamic parameters in relationship with cell culture in the shallow bubble column bioreactor. A new correlation, the turbulent energy production (TEP) of bubble-liquid two-phase flow, was established to verify the higher cell death at jetting inlet regions. Here, the effects of bubble breakup or coalescence are being neglected. The purpose of this work is to analyze the effects of H/D ratios on bubble-liquid hydrodynamics, i.e., normal and shear stress of bubble flow, energy dissipation rate of bubble phase, turbulent kinetic energy of bubble flow and two-phase turbulent energy production values, etc. The experimental observation that higher cell death at jetting inlet regions is elaborated on as well. In addition, the proposed model is based on the anisotropic characteristics instead of isotropic behaviors built in ANSYS Fluent software.

Governing Equations of the Bubble-Liquid Flows
A second-order moment bubble-liquid two-phase turbulent model based on a two-fluid continuum approach, first proposed by Zhou et al. [37,38], was improved to simulate the hydrodynamics. In this model, the bubble was considered as a dispersed phase and the liquid as a continuous phase. The bubble-bubble interactions were neglected. Single bubble motion and liquid Navier-Stokes equations under the Eulerian coordination system were established for continuity, as well as the momentum governing equations of bubble and liquid phases. Furthermore, the anisotropic behaviors of two-phase turbulent flows were fully considered. Here, we only listed the key equations and others can be referred to by the above references. They are as follows: The momentum equations of bubble and liquid phases are: The Reynolds stress transport equations of bubble and liquid phases are: where the terms on the right-hand side of Equation (2) represent diffusion term, shear production term, pressure-strain correlation term, dissipation term and phase interaction terms, respectively. These terms are provided in detail: The turbulent kinetic energy equations of the bubble and liquid phases are given: The liquid dissipation rate equation of the turbulent kinetic energy is given: where: The key problem is to close the Reynolds stress equations by a perfect anisotropic correlation to reflect anisotropic particle dispersion velocity. An interaction correlation term to successfully explain anisotropic bubble-liquid two-phase turbulence flows was first established by Zhou et al. [37]. It is: The terms on the right-hand side are closed: Thus, bubble-liquid turbulent flow can be described by the aforementioned transport equations and closure transport correlations.

Simulation Algorithm
The hybrid difference scheme to discrete the equations 1 through 13 was utilized, and the finite difference equations (FDEs) were obtained by the integration of the control volumes. The domain of calculation was first divided into a finite number of control volumes and then the differential equations were integrated over this certain control volume. A staggered grid technique was utilized to store the scalar parameters, i.e., the volume fraction, the density and the turbulent kinetic energy in the main grid points at the center of the control volumes. The tensor velocity components are marked at the control volume surface. The semi-implicit pressure linked equations algorithm, the tridiagonal marching algorithm and under-relaxation algorithm were used to solve the FDEs. The convergence criteria of 1.0 × 10 −4 for residual mass sources was set to both bubble and liquid phases.

Experiments
Two-dimensional numerical simulations were carried out in the shallow bioreactor with the height-to-diameter ratios of H/D = 2.9 (1.32 m/0.44 m) and H/D = 1.4 (0.75 m/0.58 m), respectively. Here, the swirling flow along tangential direction was neglected. The bioreactor dimensions details are provided in Figure 1. The gas entrance velocity was set to 57.2 m/s using two jetting holes with diameters of 1.0 mm in accordance with the experimental velocity value by Babosa et al. [28], corresponding to the lower superficial gas velocity of 0.06 cm/s. The initial turbulent kinetic energy was 5% of k 0 = 0.004·u jet 2 and the Reynolds stresses were specified as k' = 2/3·k 0 . The shear Reynolds stress was determined by the eddy-viscosity assumption and the outlet condition was for fully developed flow regions, in which zero gradients of all variables were taken for granted. The non-slip wall conditions were used for both gas and liquid velocities and the gas Reynolds stresses were determined by the production term, including effects of wall function on near-wall grid nodes, i.e., τ τ τ = (19) Thus, bubble-liquid turbulent flow can be described by the aforementioned transport equations and closure transport correlations.

Simulation Algorithm
The hybrid difference scheme to discrete the equations 1 through 13 was utilized, and the finite difference equations (FDEs) were obtained by the integration of the control volumes. The domain of calculation was first divided into a finite number of control volumes and then the differential equations were integrated over this certain control volume. A staggered grid technique was utilized to store the scalar parameters, i.e., the volume fraction, the density and the turbulent kinetic energy in the main grid points at the center of the control volumes. The tensor velocity components are marked at the control volume surface. The semi-implicit pressure linked equations algorithm, the tridiagonal marching algorithm and under-relaxation algorithm were used to solve the FDEs. The convergence criteria of 1.0 × 10 −4 for residual mass sources was set to both bubble and liquid phases.

Experiments
Two-dimensional numerical simulations were carried out in the shallow bioreactor with the height-to-diameter ratios of H/D = 2.9 (1.32 m/0.44 m) and H/D = 1.4 (0.75 m/0.58 m), respectively. Here, the swirling flow along tangential direction was neglected. The bioreactor dimensions details are provided in Figure 1. The gas entrance velocity was set to 57.2 m/s using two jetting holes with diameters of 1.0 mm in accordance with the experimental velocity value by Babosa et al. [28], corresponding to the lower superficial gas velocity of 0.06 cm/s. The initial turbulent kinetic energy was 5% of k0 = 0.004·ujet 2 and the Reynolds stresses were specified as k' = 2/3·k0. The shear Reynolds stress was determined by the eddy-viscosity assumption and the outlet condition was for fully developed flow regions, in which zero gradients of all variables were taken for granted. The non-slip wall conditions were used for both gas and liquid velocities and the gas Reynolds stresses were determined by the production term, including effects of wall function on near-wall grid nodes, i.e., ə ϕ/əx = 0 (ϕ = ub,ul,…). Symmetric conditions were used for both the two phases əϕl/ər = əϕb/ər = 0 at the axis as well. In-house source codes were written in the Fortran 77 language with approximately 9200 statements.
Thus, bubble-liquid turbulent flow can be described by the aforementioned transport equations and closure transport correlations.

Simulation Algorithm
The hybrid difference scheme to discrete the equations 1 through 13 was utilized, and the finite difference equations (FDEs) were obtained by the integration of the control volumes. The domain of calculation was first divided into a finite number of control volumes and then the differential equations were integrated over this certain control volume. A staggered grid technique was utilized to store the scalar parameters, i.e., the volume fraction, the density and the turbulent kinetic energy in the main grid points at the center of the control volumes. The tensor velocity components are marked at the control volume surface. The semi-implicit pressure linked equations algorithm, the tridiagonal marching algorithm and under-relaxation algorithm were used to solve the FDEs. The convergence criteria of 1.0 × 10 −4 for residual mass sources was set to both bubble and liquid phases.

Experiments
Two-dimensional numerical simulations were carried out in the shallow bioreactor with the height-to-diameter ratios of H/D = 2.9 (1.32 m/0.44 m) and H/D = 1.4 (0.75 m/0.58 m), respectively. Here, the swirling flow along tangential direction was neglected. The bioreactor dimensions details are provided in Figure 1. The gas entrance velocity was set to 57.2 m/s using two jetting holes with diameters of 1.0 mm in accordance with the experimental velocity value by Babosa et al. [28], corresponding to the lower superficial gas velocity of 0.06 cm/s. The initial turbulent kinetic energy was 5% of k0 = 0.004·ujet 2 and the Reynolds stresses were specified as k' = 2/3·k0. The shear Reynolds stress was determined by the eddy-viscosity assumption and the outlet condition was for fully developed flow regions, in which zero gradients of all variables were taken for granted. The non-slip wall conditions were used for both gas and liquid velocities and the gas Reynolds stresses were determined by the production term, including effects of wall function on near-wall grid nodes, i.e., ə ϕ/əx = 0 (ϕ = ub,ul,…). Symmetric conditions were used for both the two phases əϕl/ər = əϕb/ər = 0 at the axis as well. In-house source codes were written in the Fortran 77 language with approximately 9200 statements.
. Symmetric conditions were used for both the two phases -liquid turbulent flow can be described by the aforementioned transport equations port correlations.

ulation and Experiments gorithm
ifference scheme to discrete the equations 1 through 13 was utilized, and the finite ons (FDEs) were obtained by the integration of the control volumes. The domain of first divided into a finite number of control volumes and then the differential tegrated over this certain control volume. A staggered grid technique was utilized r parameters, i.e., the volume fraction, the density and the turbulent kinetic energy points at the center of the control volumes. The tensor velocity components are ntrol volume surface. The semi-implicit pressure linked equations algorithm, the hing algorithm and under-relaxation algorithm were used to solve the FDEs. The ria of 1.0 × 10 −4 for residual mass sources was set to both bubble and liquid phases.
ional numerical simulations were carried out in the shallow bioreactor with the er ratios of H/D = 2.9 (1.32 m/0.44 m) and H/D = 1.4 (0.75 m/0.58 m), respectively. g flow along tangential direction was neglected. The bioreactor dimensions details igure 1. The gas entrance velocity was set to 57.2 m/s using two jetting holes with mm in accordance with the experimental velocity value by Babosa et al. [28], the lower superficial gas velocity of 0.06 cm/s. The initial turbulent kinetic energy 04·ujet 2 and the Reynolds stresses were specified as k' = 2/3·k0. The shear Reynolds ined by the eddy-viscosity assumption and the outlet condition was for fully egions, in which zero gradients of all variables were taken for granted. The non-slip ere used for both gas and liquid velocities and the gas Reynolds stresses were e production term, including effects of wall function on near-wall grid nodes, i.e., ə …). Symmetric conditions were used for both the two phases əϕl/ər = əϕb/ər = 0 at the use source codes were written in the Fortran 77 language with approximately 9200 bble-liquid turbulent flow can be described by the aforementioned transport equations ansport correlations.

Simulation and Experiments Algorithm
id difference scheme to discrete the equations 1 through 13 was utilized, and the finite ations (FDEs) were obtained by the integration of the control volumes. The domain of as first divided into a finite number of control volumes and then the differential e integrated over this certain control volume. A staggered grid technique was utilized alar parameters, i.e., the volume fraction, the density and the turbulent kinetic energy rid points at the center of the control volumes. The tensor velocity components are control volume surface. The semi-implicit pressure linked equations algorithm, the arching algorithm and under-relaxation algorithm were used to solve the FDEs. The riteria of 1.0 × 10 −4 for residual mass sources was set to both bubble and liquid phases.  Figure 1. The gas entrance velocity was set to 57.2 m/s using two jetting holes with 1.0 mm in accordance with the experimental velocity value by Babosa et al. [28], to the lower superficial gas velocity of 0.06 cm/s. The initial turbulent kinetic energy 0.004·ujet 2 and the Reynolds stresses were specified as k' = 2/3·k0. The shear Reynolds termined by the eddy-viscosity assumption and the outlet condition was for fully w regions, in which zero gradients of all variables were taken for granted. The non-slip ns were used for both gas and liquid velocities and the gas Reynolds stresses were y the production term, including effects of wall function on near-wall grid nodes, i.e., ə b,ul,…). Symmetric conditions were used for both the two phases əϕl/ər = əϕb/ər = 0 at the -house source codes were written in the Fortran 77 language with approximately 9200 , bubble-liquid turbulent flow can be described by the aforementioned transport equations re transport correlations.

ical Simulation and Experiments ation Algorithm
hybrid difference scheme to discrete the equations 1 through 13 was utilized, and the finite equations (FDEs) were obtained by the integration of the control volumes. The domain of n was first divided into a finite number of control volumes and then the differential were integrated over this certain control volume. A staggered grid technique was utilized e scalar parameters, i.e., the volume fraction, the density and the turbulent kinetic energy in grid points at the center of the control volumes. The tensor velocity components are t the control volume surface. The semi-implicit pressure linked equations algorithm, the al marching algorithm and under-relaxation algorithm were used to solve the FDEs. The nce criteria of 1.0 × 10 −4 for residual mass sources was set to both bubble and liquid phases. iments -dimensional numerical simulations were carried out in the shallow bioreactor with the -diameter ratios of H/D = 2.9 (1.32 m/0.44 m) and H/D = 1.4 (0.75 m/0.58 m), respectively. swirling flow along tangential direction was neglected. The bioreactor dimensions details ded in Figure 1. The gas entrance velocity was set to 57.2 m/s using two jetting holes with of 1.0 mm in accordance with the experimental velocity value by Babosa et al. [28], ding to the lower superficial gas velocity of 0.06 cm/s. The initial turbulent kinetic energy f k0 = 0.004·ujet 2 and the Reynolds stresses were specified as k' = 2/3·k0. The shear Reynolds s determined by the eddy-viscosity assumption and the outlet condition was for fully d flow regions, in which zero gradients of all variables were taken for granted. The non-slip itions were used for both gas and liquid velocities and the gas Reynolds stresses were ed by the production term, including effects of wall function on near-wall grid nodes, i.e., ə = ub,ul,…). Symmetric conditions were used for both the two phases əϕl/ər = əϕb/ər = 0 at the ll. In-house source codes were written in the Fortran 77 language with approximately 9200 s.
hus, bubble-liquid turbulent flow can be described by the aforementioned transport equations losure transport correlations.

merical Simulation and Experiments mulation Algorithm
he hybrid difference scheme to discrete the equations 1 through 13 was utilized, and the finite ence equations (FDEs) were obtained by the integration of the control volumes. The domain of ation was first divided into a finite number of control volumes and then the differential ions were integrated over this certain control volume. A staggered grid technique was utilized re the scalar parameters, i.e., the volume fraction, the density and the turbulent kinetic energy main grid points at the center of the control volumes. The tensor velocity components are d at the control volume surface. The semi-implicit pressure linked equations algorithm, the gonal marching algorithm and under-relaxation algorithm were used to solve the FDEs. The rgence criteria of 1.0 × 10 −4 for residual mass sources was set to both bubble and liquid phases.  Figure 1. The gas entrance velocity was set to 57.2 m/s using two jetting holes with ters of 1.0 mm in accordance with the experimental velocity value by Babosa et al. [28], ponding to the lower superficial gas velocity of 0.06 cm/s. The initial turbulent kinetic energy % of k0 = 0.004·ujet 2 and the Reynolds stresses were specified as k' = 2/3·k0. The shear Reynolds was determined by the eddy-viscosity assumption and the outlet condition was for fully oped flow regions, in which zero gradients of all variables were taken for granted. The non-slip onditions were used for both gas and liquid velocities and the gas Reynolds stresses were ined by the production term, including effects of wall function on near-wall grid nodes, i.e., ə 0 (ϕ = ub,ul,…).
Symmetric conditions were used for both the two phases əϕl/ər = əϕb/ər = 0 at the s well. In-house source codes were written in the Fortran 77 language with approximately 9200 ents.
r = 0 at the axis as well. In-house source codes were written in the Fortran 77 language with approximately 9200 statements.

Results and Discussion
The key hydrodynamic parameters in association with cell culture solved by this model are the shear and normal stresses of bubble phase, the turbulent kinetic energy of bubble phase and the turbulent energy dissipation rate of bubble phase. A new definition of the energy production correlation of bubble-liquid two-phase turbulence and effects of the height-to-diameter ratios on hydrodynamics are discussed details.

Grid Independence and Validation of Proposed Model
The grid independence validation was performed based on the bubble time-averaged axial velocity (see Figure 2). Here, three kinds of grid systems, coarse (99 × 44), medium (132 × 88), and fine (264 × 176) mesh, are applied and all simulations are converged at 1.0 × 10 −4 for both bubble and liquid phases. The distributions of bubble rising velocities at the height of 0.3 m using three grid

Results and Discussion
The key hydrodynamic parameters in association with cell culture solved by this model are the shear and normal stresses of bubble phase, the turbulent kinetic energy of bubble phase and the turbulent energy dissipation rate of bubble phase. A new definition of the energy production correlation of bubble-liquid two-phase turbulence and effects of the height-to-diameter ratios on hydrodynamics are discussed details.

Grid Independence and Validation of Proposed Model
The grid independence validation was performed based on the bubble time-averaged axial velocity (see Figure 2). Here, three kinds of grid systems, coarse (99 × 44), medium (132 × 88), and fine (264 × 176) mesh, are applied and all simulations are converged at 1.0 × 10 −4 for both bubble and liquid phases. The distributions of bubble rising velocities at the height of 0.3 m using three grid systems were compared. As seen from Figure 2, all of the velocity distributions have the same trends along both axial and radial directions, especially for close agreement by medium and fine grids. Thus, the medium grid scheme was adopted for the simulation because of reasonable computation time and accepted accuracy. Figure 3 showed the validation of the proposed model by the experimental results of Lin et al. [39]. Using vertical bubble rising velocity and bubble normal stress, it indicated that they are in accordance with the experimental data. Hence, the proposed model can be used for this study with good accuracy.

Results and Discussion
The key hydrodynamic parameters in association with cell culture solved by this model are the shear and normal stresses of bubble phase, the turbulent kinetic energy of bubble phase and the turbulent energy dissipation rate of bubble phase. A new definition of the energy production correlation of bubble-liquid two-phase turbulence and effects of the height-to-diameter ratios on hydrodynamics are discussed details.

Grid Independence and Validation of Proposed Model
The grid independence validation was performed based on the bubble time-averaged axial velocity (see Figure 2). Here, three kinds of grid systems, coarse (99 × 44), medium (132 × 88), and fine (264 × 176) mesh, are applied and all simulations are converged at 1.0 × 10 −4 for both bubble and liquid phases. The distributions of bubble rising velocities at the height of 0.3 m using three grid systems were compared. As seen from Figure 2, all of the velocity distributions have the same trends along both axial and radial directions, especially for close agreement by medium and fine grids. Thus, the medium grid scheme was adopted for the simulation because of reasonable computation time and accepted accuracy. Figure 3 showed the validation of the proposed model by the experimental results of Lin et al. [39]. Using vertical bubble rising velocity and bubble normal stress, it indicated that they are in accordance with the experimental data. Hence, the proposed model can be used for this study with good accuracy.    Figure 4 and Figure 5 showed the distributions of bubble normal and shear stresses using two jetting holes (rjet/R = ±0.36) with height-to-diameter ratios of H/D = 2.9 and H/D = 1.4 at a gas entrance velocity of 57.4 m/s, respectively. Although the time-averaged velocity of bubble was relatively smaller, the higher stresses that originated from the high turbulence intensity of the liquid phase were higher because the fluid fluctuations were extraordinarily vigorous. In Figures 4a and  5a, the larger velocity gradient was observed due to stronger turbulent kinetic energy. The normal stress values gradually increased along with flow direction as well. When peaks reached the maximum values, they become a little flat due to the stable variation of stresses and velocity profiles along axial direction. Similar to those of the normal stresses variation trend, shear stresses with two   Figures 4a and 5a, the larger velocity gradient was observed due to stronger turbulent kinetic energy. The normal stress values gradually increased along with flow direction as well. When peaks reached the maximum values, they become a little flat due to the stable variation of stresses and velocity profiles along axial direction. Similar to those of the normal stresses variation trend, shear stresses with two peak values were also found in Figures 4b and 5b. When comparing their absolute values between H/D = 2.9 and H/D = 1.4 ratios, normal stresses and shear stresses had the same order of magnitude and the normal stresses were slightly larger. It indicated that the intensities of turbulent mixing and diffusion along with axial and radial directions were approximately equivalent, in which those of axial direction were slightly larger. Bubble normal stresses at axial and radial directions had larger values in the center of the reactor than those of near wall regions, which indicated that those of regions had dominantly anisotropic characteristics. Compared with the results of H/D = 2.9 to H/D = 1.4, both absolute normal and shear stress of H/D = 2.9 were greater than those of H/D = 1.4. It can be explained that the larger H/D ratio reactor easily generated more smaller-size bubbles, resulting in serious bubble breakup or coalescence, circulation and bubble fluctuations.

Turbulent Kinetic Energy of Bubble Phase
The distributions of turbulent kinetic energy of bubble phase using two jetting holes (rjet/R = ±0.36) with height-to-diameter ratios of H/D = 2.9 and H/D = 1.4 at a gas entrance velocity of 57.4 m/s are shown in Figure 6. The two peaks between center axis and reactor wall were observed and their absolute values gradually increased with the development of turbulent flow. Stronger turbulent fluctuations are always accompanied with bubble formation, movement and breakup or coalescence, which may cause the positions of peaks above jetting holes. Furthermore, according to the gas-bubble two-phase theory, non-uniform liquid velocity distributions will generate when bubbles enter the reactor through jetting holes. Larger liquid velocity gradient leads to the strong liquid fluctuations, which frequently occurred near jetting regions. The turbulent kinetic energy of the bubble phase at jetting inlet regions was weaken, which was mainly controlled by the lower bubble velocity and absence of bubble wake or bubble circulation. In comparisons with the values of H/D = 2.9, results of H/D = 1.4 were less because of the relative weak bubble and liquid fluctuations. Thus, the turbulent kinetic energy parameter is deemed inadequate to explain cell death at the jetting inlet region. The most important finding is that bubble stresses have generally smaller values near the inlet region when compared with those of fully developed flow regions. This means that normal and shear stresses increased toward the development of bubble-liquid two-phase flow, and minimum values are located near the gas entrance region. These results are consistent with experimental observations by Babosa et al. [28]. Shear stresses are contributed to by the enhancement mixing and heat and mass transfer, which are most likely to cause cell damage and death [40]. Therefore, according to the stress-associated cell damage possibility, the protentional high cell damage would occur towards fully developed flow, rather than near gas jetting regions, in terms of the simulated normal and shear stresses.

Turbulent Kinetic Energy of Bubble Phase
The distributions of turbulent kinetic energy of bubble phase using two jetting holes (rjet/R = ±0.36) with height-to-diameter ratios of H/D = 2.9 and H/D = 1.4 at a gas entrance velocity of 57.4 m/s are shown in Figure 6. The two peaks between center axis and reactor wall were observed and their absolute values gradually increased with the development of turbulent flow. Stronger turbulent fluctuations are always accompanied with bubble formation, movement and breakup or coalescence, Energies 2019, 12, 2269 9 of 16 which may cause the positions of peaks above jetting holes. Furthermore, according to the gas-bubble two-phase theory, non-uniform liquid velocity distributions will generate when bubbles enter the reactor through jetting holes. Larger liquid velocity gradient leads to the strong liquid fluctuations, which frequently occurred near jetting regions. The turbulent kinetic energy of the bubble phase at jetting inlet regions was weaken, which was mainly controlled by the lower bubble velocity and absence of bubble wake or bubble circulation. In comparisons with the values of H/D = 2.9, results of H/D = 1.4 were less because of the relative weak bubble and liquid fluctuations. Thus, the turbulent kinetic energy parameter is deemed inadequate to explain cell death at the jetting inlet region.

Turbulent Kinetic Energy of Bubble Phase
The distributions of turbulent kinetic energy of bubble phase using two jetting holes (rjet/R = ±0.36) with height-to-diameter ratios of H/D = 2.9 and H/D = 1.4 at a gas entrance velocity of 57.4 m/s are shown in Figure 6. The two peaks between center axis and reactor wall were observed and their absolute values gradually increased with the development of turbulent flow. Stronger turbulent fluctuations are always accompanied with bubble formation, movement and breakup or coalescence, which may cause the positions of peaks above jetting holes. Furthermore, according to the gas-bubble two-phase theory, non-uniform liquid velocity distributions will generate when bubbles enter the reactor through jetting holes. Larger liquid velocity gradient leads to the strong liquid fluctuations, which frequently occurred near jetting regions. The turbulent kinetic energy of the bubble phase at jetting inlet regions was weaken, which was mainly controlled by the lower bubble velocity and absence of bubble wake or bubble circulation. In comparisons with the values of H/D = 2.9, results of H/D = 1.4 were less because of the relative weak bubble and liquid fluctuations. Thus, the turbulent kinetic energy parameter is deemed inadequate to explain cell death at the jetting inlet region.

Turbulent Energy Dissipation Rate of Bubble Phase
The definition of the turbulent energy dissipation rate of bubble phase is: The distributions of turbulent kinetic energy of the bubble phase with height-to-diameter ratios of H/D = 2.9 and H/D = 1.4 at a gas entrance velocity of 57.4 m/s are shown in Figure 7. We can see that the ε began to increase toward the development of turbulent flows above the jetting holes. The largest energy dissipation rate of bubble flow located at fully developed regions can be explained by Kolmogorov's energy cascade theory, that energy is dissipated by fluctuating viscous stresses in resisting deformation of fluid at fluctuating strain rates. Hence, a dissipation rate is primarily

Turbulent Energy Dissipation Rate of Bubble Phase
The definition of the turbulent energy dissipation rate of bubble phase is: Energies 2019, 12, 2269 10 of 16 The distributions of turbulent kinetic energy of the bubble phase with height-to-diameter ratios of H/D = 2.9 and H/D = 1.4 at a gas entrance velocity of 57.4 m/s are shown in Figure 7. We can see that the ε began to increase toward the development of turbulent flows above the jetting holes. The largest energy dissipation rate of bubble flow located at fully developed regions can be explained by Kolmogorov's energy cascade theory, that energy is dissipated by fluctuating viscous stresses in resisting deformation of fluid at fluctuating strain rates. Hence, a dissipation rate is primarily determined because of rapid viscous scale dissipation, and energy is set down by the non-linear process of scale-to-scale energy. Considering the effects of the larger H/D ratio, a larger energy dissipation rate was observed in comparison with the small H/D ratio because much more smaller bubbles originating from bubble breakup or coalescence caused severe bubble interactions and circulation and brought higher energy dissipation consumptions. In the meantime, the larger energy dissipation rate at fully developed regions can be elaborated by this explanation. The bubble energy dissipation rate did not exhibit large values at the jetting inlet region, and cells also located in the liquid phase were coupled with the low liquid phase energy dissipation rate. They indicated that the energy dissipation rate parameter may not be a satisfactory hydrodynamic parameter to reveal cell damage at jetting inlet regions. In summary, the aforementioned hydrodynamic parameters, i.e., the normal and shear stresses, the turbulent kinetic energy and the turbulent energy dissipation rate, indicated that all of them have smaller values at the gas inlet region when compared to those of the developed regions. Therefore, they may not be responsible for the higher cell death rate observed at the jetting inlet region. Hence, an appropriate parameter that can theoretically explain such a phenomenon is In summary, the aforementioned hydrodynamic parameters, i.e., the normal and shear stresses, the turbulent kinetic energy and the turbulent energy dissipation rate, indicated that all of them have smaller values at the gas inlet region when compared to those of the developed regions. Therefore, they may not be responsible for the higher cell death rate observed at the jetting inlet region. Hence, an appropriate parameter that can theoretically explain such a phenomenon is required.

Turbulent Energy Production of Bubble-Liquid Two-Phase Flow
So far, a reasonable hydrodynamic parameter that successfully discloses the higher cell death rate near the jetting inlet region has not been obtained. Here, we define a correlation of the turbulent energy production (TEP) of bubble-liquid two-phase flow to perform a theoretical investigation. It is: The turbulent energy production correlation is a function of the Reynolds stress and the mean velocity gradient, representing the energy exchange rate between mean and fluctuating motions. It is quietly different from the aforementioned energy dissipation that presented by the form of viscous heat loss. Figure 8 shows the bubble-liquid two-phase turbulent production term contour near the inlet region and fully developed regions using two jetting holes (rjet/R = ±0.36) with H/D = 2.9 and H/D = 1.4 at a gas entrance velocity is 57.4 m/s, respectively. As we can see, the values near the jetting inlet regions are significantly larger than those of the fully developed regions due to the larger energy transfer rate in that region coming from the gas injection. It is quietly different from the energy dissipation rate, characterized by the viscous heat loss. These results were satisfied well with the experimental results [28]. A unimodal peak was observed between two jet holes, which may have originated from the bubbles' synergistic effects by the neighbor jetting holes. Meanwhile, bimodal distributions near the jetting region, and flatter profiles at the developed flow regime, were found. It can be explained that when stronger bubble interaction effects happened, they occurred among multi-holes near the jetting inlet regions, and gradually weakened with the development of flows. The H/D ratio values near the jetting inlet regions were greater than those of the developed regions since the bubbles attached to the sparger surface, which installed at the reactor's bottom section, had larger bubble sizes as generation during the first stage. They then begin to decrease with the bubble flow toward the top of the reactor together with bubble swing, downward flow, circulation, deformation and bubble plume fluctuations.
Compared with H/D = 1.4, those values of H/D = 2.9 are greater (see Figure 8). As a matter of fact, the bubble behaviors near the jetting inlet regions belong to the sparger configuration impact regions on bubble-liquid two-phase flow regime. The bubbles generated at sparger surfaces began to grow up toward the developed flow gradually, accompanied with coalescence, breakup, swing and circulation fluctuations. It enhanced the energy dissipation and weakened the fluctuations along with the developing flows. So, these values at a far distance from the sparger position should be smaller than those near the jetting inlet regions. Furthermore, the smaller H/D ratio under a fixed volume condition, representing a larger diameter, was beneficial for the growth of bubble expansion, the mitigation of bubble rupture and bubble-liquid fluctuation intensity. Thus, the Favre-averaged Reynolds stresses caused by fluctuations were reduced, which is in favor of the smaller values generation.
According to Figures 8 and 9, the maximum energy production values were approximately 40.0 W/m 3 , much less than the approximate 3000 W/m 3 simulated by Liu et al. [27] using Babosa's experiment reactor dimension (see Figure 10). This extreme narrow-shape bioreactor, with an H/D ratio of 20.1, has bigger criteria threshold value than that of H/D = 5 [9,11], which gives rise to the largest values near the jetting inlet regions resulting in the deterioration of cell culture. circulation fluctuations. It enhanced the energy dissipation and weakened the fluctuations along with the developing flows. So, these values at a far distance from the sparger position should be smaller than those near the jetting inlet regions. Furthermore, the smaller H/D ratio under a fixed volume condition, representing a larger diameter, was beneficial for the growth of bubble expansion, the mitigation of bubble rupture and bubble-liquid fluctuation intensity. Thus, the Favre-averaged Reynolds stresses caused by fluctuations were reduced, which is in favor of the smaller values generation. According to Figures 8 and 9, the maximum energy production values were approximately 40.0 W/m 3 , much less than the approximate 3000 W/m 3 simulated by Liu et al. [27] using Babosa's experiment reactor dimension (see Figure 10). This extreme narrow-shape bioreactor, with an H/D ratio of 20.1, has bigger criteria threshold value than that of H/D = 5 [9,11], which gives rise to the largest values near the jetting inlet regions resulting in the deterioration of cell culture.

1.
The proposed second-order moment bubble-liquid two-phase turbulence model can well simulate the hydrodynamic parameters with good accuracy.

2.
The velocities and their fluctuations, and the normal and shear Reynolds stresses of bubble and liquid phases are distinctively anisotropic.

3.
The established turbulent energy production successfully verifies the higher cell damage at the gas jetting inlet region. 4.
The shallow shape bioreactor is more favorable for cell living than the long and narrow-shaped reactor.
Author Contributions: All authors equally contributed to the data analysis and the simulation, the results analysis, the writing, and review