Development of modeling and simulation of bubble‐liquid hydrodynamics in bubble column

In the present work, the investigation on strategies to alleviate the cell death in terms of bubble‐liquid two‐phase hydrodynamics in bubble column was performed using the modeling and numerical simulation approaches. A multifluid bubble‐liquid turbulent model based on the second‐order moment algorithm was developed to predict turbulent hydrodynamics. In this model, the anisotropic characteristics of Reynolds stresses and the interactions between bubble and liquid phases are considered by means of the defined bubble‐liquid two‐phase fluctuation velocities correlation. Hydrodynamic parameters such as turbulent energy dissipation rate, the turbulent kinetic energy, the distributions of bubble normal and shear stresses, and the operation conditions of gas entrance velocity and height to diameter of column are accurately simulated. A proposed correlation named turbulent energy production rate is revealed successfully the higher cell death rate at the gas inlet regions that cannot be verified by the rest of the two‐phase hydrodynamics. Superficial gas velocity is the predominant influence on this value distribution as smaller the ratio of height to diameter of bubble column, which reduced by two order of magnitude under decreased by one level of gas velocity. Near inlet region with large height‐to‐diameter ratio, this value is always larger than those of small one.

stress-sensitive parameter that induced by bubble behaviors can be considered as one of the key causes for cell damage. Jöbses et al 5 categorized the three kinds of bubble behaviors in bubble column reactor as pertaining to potential cell damage, that is, bubble formation at near gas inlet region, bubble rising region at developed flow regions, and bubble rupture at the top liquid interface. Most researcher believed that bubble breakup near liquid interface is the mainly contribution to cell death. Babosa et al 6,7 found that bubble formation at the sparger was the main event leading to cell death rather than bubble rising and breakup during Dunaliella tertiolecta and Dunaliella salina cultivated in bench-scale bubble column. It consists of two stages: In the first stage, bubbles will expand with sparge surface while remaining being attached to nozzle. In the second stage, the buoyance force incurs the bubble to detach from nozzle. When cells are in contact directly with bubbles or near sparger, the shear force caused by rapidly expanding bubbles and the two-phase turbulent flow caused by liquid rushing into cavity formed by bubble detachment can be potential source for cell damage. 8 Camacho et al 9 reported that microalgae Phaeodactylum tricornutum cells are sensitive to hydrodynamic stress, which is caused by the interaction between bubbles, and they attributed the cause of cell damage factors to the breakup of smaller bubbles on liquid surface. Therefore, the maximum hydrodynamic stress in complex turbulence plays very important roles to characterize the threshold values for cell death. However, it is very difficult to obtain this threshold value due to unclearly mechanism understanding of multiphase turbulent flow.
Early studies have usually focused on the laminar flow mode system, and the bubble behaviors were not be considered. [10][11][12][13][14][15] Hu et al 16 found that threshold values for mammalian cells growing in suspension cover a wide range of from 8 to 450 Pa and cannot obtain a precise value. Nienow et al 17 indicted that if hydrodynamic stress is lower than critical value, it makes no effect on cell growth; however, the negative response will result in the slower growth, change of metabolism, etc, when it exceeds over threshold value. Sieck et al 18 revealed that the oscillating hydrodynamic stress will aggravate effect on cell damage in larger-scale bioreactor. Another effective analysis method, the so-called energy dissipation rate and eddy-length-scale hydrodynamic parameters, can be used to quantify the cell damage that caused by bubble burst. Cherry et al 19 and Croughan et al 20 disclosed that cell injury is most likely events if the smallest energy containing eddy has the same size or smaller than the cells. Christi 1 proposed the Kolmogoroff microscale-eddy length correlation and also argued that rupturing bubbles are enough to damage cells. Ma et al 10 indicated that cell damage may occur when the local energy dissipation rate (EDR) coming from small bubble breakup was >10 000 kW/m 3 when the eddy-length scale reaches typical cell diameter of 10-20 μm. Particularly, small size bubbles (<2 mm diameter) are more dangerous to kill cell during bubble rupture in comparison with larger bubble (~10 mm diameter). Mollet et al 11,21 suggested that the primary cause for cell damage was necrotic under higher EDR levels based on microfluidic device in which cells are subjected to well-defined hydrodynamic forces, as well as the apoptosis was also revealed at lower and sublysis EDR levels. Hu et al 16 summarized all the key research findings on EDR and its implications on cell growth. Generally, positive relations between EDR and cell damage were found. Liu et al 22 proposed a new hydrodynamic correlation that related to cell culture that he considered the stress was induced by the two-phase turbulent energy production.
The reactor design, particularly the ratios of height to diameter take great effect on cell living. Suzuki et al 23 found a linear correlation between specific death rate and inverse of culture height. However, the inverse trend that an increase in culture height resulted in an increased death rate is observed by Ref. 24 The larger height will lead to capture more cells by the rising bubbles that carried to liquid surface, where cells die as the bubbles rupture was related. Moreover, hydrodynamic conditions will not promote cell damage only at the laboratory scale, which could not be extended to the pilot scale when using similar spargers and superficial gas velocities. 9 The effects of height-to-diameter (H/D) ratio are found to be important due to flow transition from homogeneous to heterogeneous regime and vice versa. 9,25-27 Throat et al 28 investigated the effect of combined sparger design and H/D ratio on gas holdup. They found that gas holdup has a maximum at H/D ratio of 1 and decreases by 15%-20% as the H/D ratio increased to a value between 4 and 5 for a multipoint sparger in an air-water system. However, for single point sparger, the reverse trend was found by Lau et al, 29 which used three types of gas distributors (single nozzle, perforated plate, and porous plate) to study the effect of H/D ratios between 2 and 7.2 on hydrodynamics. An increase in H/D ratio resulted in a decrease in the overall gas holdup, with a diminishing effect at values higher than H/D ratio of 4. Veera et al 30 found that gas holdup increased by about 50%-100% as H/D increase from 1 to 5 for a single point sparger and it reaches the maximum at very low H/D and decreases by 15%-20% as H/D increases up a value of 5 for a multipoint sparger. Dhotre et al 31 found that centerline gas holdup increased and wall gas holdup decreased with distance from bottom increase for multipoint sparger and gas holdup profiles are very steep at the axial location of H/D = 0.259. If the homogeneous regime is ensured that under the condition that the diameter is <150 mm, superficial gas velocity is <100 mm/s, and the H/D is larger than 5, the influence becomes negligible. As mentioned above, the knowledge or criteria regarding minimized cell damage as scaling up has still unclearly. Even though the experimental method has been successfully conducted as a primary research way, they still being limited on the laboratory scale and the prediction parameters that needed to be constructed are quietly confined. Using theoretical approach coupled with accurate correlations and validated by experiments, pilot-plant and industrial-scale bioreactors can be potentially constructed.
The computational fluid dynamics became a promising approach that uses a numerical simulation to predict the single-phase/multiphase laminar and turbulent flow conditions, that is, gas flow, bubble-liquid two-phase flows, gas-particle two-phase flows, and gas-liquid-particle three-phase flows. Mollet et al 21 conducted a simulation using ANSYS Fluent to predict the EDR distributions in a 1.6-L bioreactor with a standard Rushton turbine impeller at two impeller speeds (73 and 195 rpm). They found that the local energy dissipation rate scale is on the order of the 10 7 W/m 3 ; this resulted in an immediate, catastrophic rupture of suspended animal cells that the averaged energy dissipation rate is on the order of the 10 3 . Rajaram et al 32 performed a modeling and simulation using a standard k-ε two-equation turbulence model built-in ANSYS Fluent to predict the distributions of average energy dissipation rate, maximum energy dissipation rate, average shear rate, and average normal stress of stirred tank bioreactors. They found that tank diameter, impeller diameter, width, angle, and blade type effects on the extent of deactivations were dramatically different. There was also a strong correlation with the average turbulent normal stress and the mass transfer coefficient. Zhang et al 33 proposed multiple size groups (MUSIG) model based on the population balance model and the Eulerian-Eulerian two-fluid approach, which consider the effects of bubble size distributions on hydrodynamics. The reactor parameters in their study included a stirred 15.0-L bioreactor with three six-bladed Rushton impellers and four equally spaced baffles. The mass transfer varied with the smallest sized bubbles around 2 mm in diameter, and highest mass transfer located near the impeller region. Chris et al 34 considered that turbulent eddies are much larger or smaller than the biological cell diameter, which would not inhibit or damage cell viability. However, once eddies are similar in size to the cell diameter, they become harmful. Based on initial and refined simulations done with a standard k-ε, two-equation turbulence model, and the sevenequation Reynolds stress model using ANSYS Fluent, the energy dissipation rates and the Kolmogorov length-scale distributions for the aerated stirred tanks are evaluated. The results showed that the mean turbulent eddies in all scales are substantially larger than the average cell diameter, indicating that the cell death due to mechanical agitation exists at all scales.
In this work, a second-order moment bubble-liquid twophase model, considering completely anisotropic characteristics of bubble dispersion, was proposed to simulate the hydrodynamic parameters in bubble column details. The strategies to alleviate cell damage, that is, lower gas entrance velocity and different height-to-diameter ratios of bioreactor, were investigated, particularly verification for higher cell death near gas inlet region by proposed the conception of turbulent energy production rate.

| GOVERNING EQUATIONS OF BUBBLE-LIQUID T WO-PHASE TURBULENT FLOW
In this study, we have developed a second-order moment (SOM) bubble-liquid two-phase turbulent model based on a two-fluid continuum approach, firstly proposed by Zhou et al, 35,36 in which gas bubble is treated as dispersed phase and liquid is as continuous phase. Continuity and momentum governing equations of liquid-phase and bubble-phase under Eulerian coordinates are established based on the single-bubble motion equation and the liquid-phase Navier-Stokes equations. The anisotropic behaviors in two-phase turbulent flows are considered completely in this model, which is better than the traditional isotropic hypothesis for bubble dispersions. Both gas and liquid volume fractions vary to a great extent for bubble-liquid turbulent flows. In order to reduce the number of closure equations, mass weighting and time averaging are applied to the governing equations. Generalized forms of continuum and momentum equations, as well as Reynolds stress transport equations of bubble and liquid phase, and interaction closure transport equation are listed below: The continuity equations of bubble and liquid phases are The momentum equations of bubble and liquid phases are The two phases governing equations were derived from the Reynolds stress equations for the single-phase flows. The first is obtained by multiplying the instantaneous volume averaged momentum equations for two velocity components in the i and j directions by the velocity components at the j and i directions. The second is acquired by multiplying the averaged momentum equations for the two velocity components in the i and j directions by the averaged velocity components at the j and i directions. Subtracting one equation from the other, the Reynolds stress transport equation for two-phase flows can be obtained as where the terms on the right-hand side of Equation (3) represent diffusion term, shear production term, pressure-strain correlation term, phase interaction term, and dissipation terms, respectively. The detailed forms of these terms are given as below: The turbulent kinetic energy equations of the bubble and liquid phases are The liquid dissipation rate equation of the turbulent kinetic energy is where As a matter of fact, the current correlation is a typically isotropic velocity correlation that derived from nondimensional analysis, which is not satisfied match with experimental results. Regarding bubble-liquid two-phase turbulent flow, we developed a new correlation to consider successfully the anisotropic characteristics, firstly established by. 35 The closure transport is listed as Equation (14) The terms of right-hand side are closed by the following set of correlations.
In summary, all the governing equations above have been closed including a set of mass and momentum conservation equations, their Reynolds stress, and closure transport equations.

| 1 Algorithm
The above-governing equations were discretized using the hybrid difference scheme to a set of the finite difference equations. The calculation domains were firstly divided into the control volumes and then were integrated over this certain control volume. The volume fraction, the density, and the turbulent kinetic energy are stored the main grid points located the center of control volumes, the velocity components are identified at the control volume surface by means of the staggered grid method. Finite difference equations are solved by the combination of the semi-implicit pressure linked algorithm that used to correct the p-v correlation, the tri-diagonal marching algorithm line-by-line iteration, and under-relaxation algorithm. The convergence values of 1.0 × 10 −4 are set to the residual mass sources of bubble and liquid phases. The in-house written by Fortran 95 with approximately 11 500 statements.

| 2 Reactor dimensions and boundary conditions
The 2D bubble column consists of the diameter of 14.9 cm and the height of 160 cm. Two ratios of height to diameter (H/D) are set to 6.6 and 3.3. The detailed simulation parameters for operation condition and reactor dimensions are given in Figure 1. As for boundary conditions, two kinds of superficial gas velocities are the 1.0 and 0.5 cm/s, initial turbulent kinetic energy was specified as 5% of the square of velocity k 0 = 0.005 * u jet 2 , and Reynolds stresses were specified as 2/3 of the turbulent kinetic energy k′ = 2/3 * k 0 . The fully developed flow conditions of two phases were set to as outlet conditions, which zero gradients of all variables were taken. Nonslip wall condition was used for both gas and liquid velocity, and the gas Reynolds stresses were determined by the production term including the effect of wall function for nearwall grid nodes, that is, əϕ/əx = 0 (ϕ = u b , u l , …). At the axial direction, symmetric conditions are adopted for both the two phases əϕ l /ər = əϕ b /ər = 0.

| RESULTS AND DISCUSSION
Bubble movement begun to occur once the jetting gas enters water and then give rise to the upflow toward surface as a result of buoyance and drag forces. All solutions are satisfied with convergence, and the hydrodynamic parameters are solved by the proposed multifluid model and source code.

| Grid independence test
The independence test of grid size resolution based on the sensitivity of bubble time-averaged axial velocities to spatial is conducted (see Figure2). Here, the coarse system of 200 × 125, the medium system of 400 × 250, and fine system of 600 × 360 mesh grids with axial by radial are used. We can see that the velocity distributions have similar trends along both axial and radial directions, especially for medium and fine grids. Therefore, the medium grid size scheme was utilized due to economic computation time and reliable accuracy.

| Simulations of bubble-liquid hydrodynamic parameters
The proposed multifluid turbulent model to simulate the key hydrodynamic parameters associated with cell culture was conducted. Two kinds of lower superficial gas velocity of 1.0 and 0.5 cm/s and two height-to-diameter ratios of H/D = 6.6 and H/D = 3.3 were adopted in the simulation cases. Although the time-averaged velocity of the bubble is relatively small and their values generally <20 cm/s (see Figure 3), the higher stresses that caused by strong turbulence intensity relative to the liquid phase are observed. Three distinct peaks of bubble axial normal stresses between the central axis and the reactor wall above the gas entrance region are found, in which the larger velocity gradients were generated due to larger turbulent kinetic energy. With increasing along with axial position toward the fully developed region, normal stress values gradually increased, and when it reaches up to peaks, they become slightly flatter due to much more stable changes of stresses and velocity profiles in axial location (see Figure 4A,B). The same trend of shear stress distributions with two peak values can also be indicated in Figure  4C,D. The larger superficial gas velocities have corresponded with bigger stress values than those of smaller superficial gas velocities due to larger input of energy and turbulent intensity. Another observation is that bubble stresses are generally smaller at near inlet region when compared with the fully developed flow regions. These means that normal and shear stresses increased with the development of bubble-liquid twophase flows. As shown in Figure 4C,D, bubble shear stresses are always smaller at jetting inlet that caused by the large velocity gradients near the jetting flow regions wall and tend to decrease with increasing radial positions. Shear stress is a main reason for enhancing mixing and heat and mass transfer, which is also likely to cause cell damage and death, by Ref. 4 Compared to normal stresses, the absolute values of shear stresses are smaller and similar trends at other directions as well. Because of higher viscous dissipations caused by wall boundary effect, both normal and shear stresses of H/D = 3.3 are smaller than those of H/D = 6.6, even they have similar trend along axial development as well. Therefore, high cell damage would occur toward fully development flow instead of at near jetting inlet region.

| Distributions of second stress tensor of bubble phase
Miguel et al 37 indicated that stress alone either axial direction or radial direction is not adequate to give threaten to cell damage under complex flow conditions. He proposed an invariant second stress tensor as a hydrodynamic parameter to characterize and predict cell damage. It is The problem of this definition is that it neglects the correlations between normal stress and shear stress.
Here, we improved it as below The comparisons of ratio and stress tensor considering and without considering this correlations between normal stress and shear stress are shown in Figure 5A,B. We can see that the both their values at the gas inlet are smaller, as well as those values considering correlation are smaller than those of without considering it. Figure 6  It can be seen that tensors are smaller at the gas inlet region when compared to the moving developed flow region. Under smaller superficial gas velocity and smaller H/D values, these values are always smaller and both have similar trend regarding normal and shear distributions. Therefore, this definition also fails to explain the higher cell death rate at the gas entrance region.

| Distributions of turbulent energy dissipation rate of bubble phase
The turbulent energy dissipation rate of bubble phases is defined as Figure 7 shows the distributions of bubble turbulent kinetic dissipation rate under various H/D ratios and (21) Since energy dissipation rate is controlled by the energy transfer rate coming from larger energetic scales, it is always >0. Smaller superficial gas velocities correspond to smaller energy dissipation, similar to those of stress distributions trends as well. However, larger energy dissipation rate at H/D = 3.3 is a bit larger than that of H/D = 6.0. It can be explained that the high dissipation rate due at smaller height-to-diameter ratio resulted in much more small size bubbles originated from bubble breakup or coalescence, which causes the severe bubble interactions, leading to higher energy dissipation. Bubble-phase energy dissipation rates have still smaller at near inlet region. In addition, cells located in this region were also were accompanied with the lower liquid-phase energy dissipation rate. Both lower energy dissipation rate of bubble and liquid phases demonstrated that it may not be an ideal way to explain higher cell damage.

| Distributions of turbulent kinetic energy of bubble and liquid phases
The distributions of turbulent kinetic energy of both liquid and bubble phases are shown in Figures 8 and 9. As can be seen from Figure 8A,B, two peaks located between center axis and wall are observed and they gradually increase with superficial gas velocities. Based on the gasbubble two-phase turbulence theory, nonuniform liquid velocity distribution was generated as bubbles enter the reactor through jet holes. Larger liquid velocity gradient led to the generation of liquid turbulence at near jetting regions. Meanwhile, liquid turbulent fluctuations are also induced by bubble wake formation, bubble swing, and bubble downward circulation flows. The smaller values at inlet regions are a result of lower bubble velocity and the absence of bubble wakes, and the larger values at near-wall region are caused by wall region effects, stronger turbulent mixing and diffusions induced stronger turbulence. As a matter of the fact that larger superficial gas velocity has the bigger kinetic energy that came from the larger kinetic energy mean flows and interactions among bubbles such as coalescence of smaller bubbles, formation of big bubbles, and breakup of larger bubbles. In Figure 9A,B, three peaks are obtained between the center axis and the wall at near gas inlet and it gradually increased with the development of turbulent flow and superficial gas velocity. These phenomena can be attributed to the effects of lower bubble velocity, absence of bubble wakes, bubble swing, and downward circulation flow. The stronger turbulent fluctuations are always found that they were accompanied with bubble formation, movement, and breakup or coalescence, resulting in the peaks located above jetting holes. Furthermore, we can see that the values of larger H/D are always smaller than those of smaller H/D because that the high proportion of small bubble sizes was easier to generate at a higher height of H/D = 3.3 due to bubble breakage and coalescence mechanism, which gave more contributions to the vigorous turbulent fluctuation and intense kinetic energy transfer. Compared with bubble phase, the liquid turbulent kinetic energy is slightly smaller, and they showed similar trends such as smaller values at near inlet region. This can be caused by the effects of lower bubble and liquid velocities; the absent bubble wakes at near inlet regions. Hence, turbulent kinetic energy is inadequate to explain the higher cell death at the jetting inlet. A summary of the effects of above-mentioned hydrodynamic parameters is that they have smaller values at the gas inlet region when compared to the developed region. Therefore, they may not be the reason for interpretation of the higher cell death rate at the jetting inlet region. Hence, it is necessary to build up a reasonable correlation or parameter to explore this phenomenon.

| Distributions of turbulent energy production correlation of bubble phase
Mainly, investigations regarding cell damage are attributed to the shear stress production, suggested by Charmers et al. 2 In this study, we have proposed a turbulent energy production correlation (TEPC) based on the two-phase theory to elaborate the higher death rate at near gas inlet region. Essentially, it is a derivation of bubble shear production. This correlation was defined by Figure 10 indicates the contour distributions of this values (W/m 3 ) under H/D = 6.6 and u sg = 1.0 cm/s, respectively. We can see that those at near inlet region are significantly larger than fully developed region regions, which may be represented that large energy transfer rate in that region resulting from gas injection. This definition is quietly different than the energy dissipation rate, which is characterized by the viscous heat loss. Another notable phenomenon is that the maximum peaks with unimodal distributions were observed at between the two jet holes. This may result from the synergistic effect of the bubbles generated from adjacent jet holes. The larger values are founded due to an increase in gas entrance velocity that carried more energy input and intensified process of bubble breakup, deformation, and coalescence. In the meantime, those at near inlet region under smaller H/D are smaller than larger H/D. It can be explained that the larger wall boundary effects caused by bigger H/D ratio had generated the more energy dissipations, resulting in the decrease in the energy transfer. However, there are almost no effects by H/D under developing flow regions. The reasons are that bubbles attached sparger surface at reactor bottom have the larger bubble size and begun to decrease along with the bubble flow toward reactor top. More bubbles begun to breakup to produce the much more energy dissipation rates. Hence, the kinetic energy transferred from turbulent eddy was reduced to lead to a decrease in TEPC. The effects of H/D ratio caused by the wall viscous dissipation may be neglected. It can be seen from Figure 12 that the maximum TEPC values in this simulation are approximately 25.0 W/m 3 , which is much less than reported predicted range, approximately 3500-4000 W/m 3 by Babosa et al. 6 This experimental conditions are the diameter of D = 3.58 cm and height H = 72.5 cm with H/D = 20.1 and u sg = 3.9-7.8 cm/s. It is obvious that this kind of extremely narrow shape of reactor is easier to generate larger TEPC values at near inlet regions, and its effect on TEPC production is much larger than that of increasing superficial gas velocity, which resulted in the exacerbation of cell death rate at near inlet region. Therefore, smaller H/D ratio and smaller superficial gas velocity are contributed to the alleviation of cell death. But, due to the limitation of production capacity that requires higher gas entrance velocity for commercial application, the reasonable selections for H/D ration and superficial gas velocity need to be further optimized. 2. The anisotropic characteristics of bubble Reynolds stresses and interaction between bubble and liquid were indicated. 3. Higher cell damage at near gas jetting was verified by proposed turbulent energy production correlation reasonable. 4. More critical transition thresholds to indicate the cell damage under different H/D ratios will be further investigated in next work.

ACKNOWLEDGMENT
We sincerely appreciate the financial support of the National