Numerical Simulation of Flow in Erlenmeyer Shaken Flask

1. Introduction

By far most of all biotechnological experiments are carried out in shaken bioreactors [1,2]. Every laboratory bioprocess development developing in early stages are relied on parallel thousands of experiments in shaken flasks to determine optimal medium composition or to find an suitable microbial strain due to the great experimental simplicity of the apparatus. Furthermore, it can also helpful for decisive and orienteering decisions on experimental conditions. However, shaken flask experiments could only provide phenomenal conditions such as rotary speed which reflects mixing and oxygen requirement degree in aeration process, it could not quantify important engineering parameters like the volumetric power consumption [3], the oxygen transfer capacity [4-6] or the hydro-mechanical stress [7,8] which would be more crucial for process scale-up. Thus such parameters have to be determined through empirical or semi-empirical equations or depended on pilot experiments. This facts is doubtfully wakened the reliability of shaken results and also prolong the period of bioprocess through laboratory to industrial process. To gain deeper understanding of the afore mentioned mechanisms on a theoretical basis, the geometry, i.e. the contour and spatial distribution of the rotating liquid mass moving inside a shaken Erlenmeyer flask is of crucial importance. The liquid distribution gives important information about the momentum transfer area, which is the contact area between the liquid mass and the flask inner wall, and the mass transfer area, which is the surface exposed to the surrounding air, including the film on the flask wall [4]. B¨uchs et.al [9] have setup the liquid distribution model flow characterization of liquid in shaken flask to calculated the liquid distribution, and validated with photography. However their calculation based simply liquid shape did not consider the liquid surface bend or sunk during rotation, thus the calculated maximum height of liquid approached has a little difference with experimental results, and also this model could not calculate the gas-liquid interface which may important for oxygen transfer.

Computational fluid dynamics (CFD) is a novel method to investigate the flow behavior with low cost, independent on container geometry. It can also provide more details which could not be obtained by experiments. Many commercial CFD softwares such as Phoenics, Fluent, CFX, StarCD have presented excellent success in simulation for both process and apparatus.

This work aims at proving dynamic fluid kinetic model with dynamic mesh of shaken flask, to calculate the liquid distribution, contact area between the liquid mass and the flask wall and energy dissipation etc. in shake flasks to provide the information required to investigate momentum and gas/liquid mass transfer and volumetric power consumption. All the work was carried out with Fluent 6.2 in this work.

2. Model formulation

2.1. Basic concepts of shaken movement

Most commonly shake flasks are agitated by orbital shaking machines. Figure. 1 illustrates the physic description of this shaking motion. The shake flask performs a circular translatoric movement with the radius equal to half the shaking diameter keeping its orientation relative to the surrounding. This movement is synthesized by a superposition of two individual movements. The first movement is the angular velocity ω ₁ of the circular translatoric movement with the shaking radius around the center of the shaking motion (shaft 1), which is related to as the shaking frequency of the shaking machine, N.

$${\omega }_1=2\pi N$$ E1

Figure 1.

Theoretical partition of shake flask movement: superposition of circular translatoric and rotational movement.

The second movement is the rotation counteracted the first rotation to keep the shake flask’s spatial orientation around the flask center (shaft 2). Therefore: ω ₂ = −ω ₁, driven by the centrifugal field of the first rotation (−ω ₁).Thus, in shaken flask, the liquid will move in ellipses in caterian coordinate system, as described as equation 1 and 2.

$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$$ E2

$$\begin{array}{l}{U}_X=ab\omega \times \frac{\Lambda \cdot \cos (\omega t)-\Pi \cdot \sin (\omega t)}{{\Lambda }^2}\\ U_Y=ab\omega \times \frac{\Lambda \cdot \sin (\omega t)-\Pi \cdot \cos (\omega t)}{{\Lambda }^2}\end{array}$$ E3

Where, ω is the angular velocity, t is the rotation time. a and b are radius in shaft 1 and shaft 2 axial respectively. Parameters Λ and Π are defined in equation (3).

$$\begin{array}{l}\Lambda =\sqrt{a^2{\cos }^2(\omega t)+b^2{\sin }^2(\omega t)}\\ \Pi =\frac{b^2-a^2}{2\Lambda }\sin (2\omega t)\end{array}$$ E4

For most Erlenmeyer shaken flask is designed as circular radius, so a equals b, the movement equation can be simplified as equation (4) and (5).

$$x^2+y^2=a^2$$ E5

$$\begin{array}{l}Ux=a\omega \cos (\omega t)\\ Uy=a\omega \sin (\omega t)\end{array}$$ E6

2.2. Model for flow dynamics

In most cases, there are two phases of both gas and liquid in shaken flask, and there is little entanglement between both of them, gas phase and liquid phase are always separated with little exchange. As the results, Volume of Fluid Model of Eularian model is more practical for the description of motion in shaken flask.

Continua equations:

In VOF model, the interface between phases is determined by the volume void of one or more phase in the multiphase system. The continua equation of phase q is as follow.

$$\frac{1}{{\rho }_q}\left[\frac{\partial }{\partial t}(a_q{\rho }_q)+\nabla \cdot (a_q{\rho }_q{\overrightarrow{\upsilon }}_q)\right]=S_a{}_q+{\displaystyle \sum _{p=1}^n(\stackrel{.}{{\dot{m}}_{pq}-{\dot{m}}_{qp}}})$$ E7

In which, ρ _q and α _q presents density and volume void for q phase respectively. p presents all phases excluding q phase in calculation region, and S _αq is the resource item. m _pq states the exchange from p phase and q phase. m _qp states the exchange from q phase and p phase. In this case, if all the phase exchanges are neglected, then equation (6) can be rewritten as equation (7).

$$\frac{1}{{\rho }_q}\left[\frac{\partial }{\partial t}(a_q{\rho }_q)+\nabla \cdot (a_q{\rho }_q{\overrightarrow{\upsilon }}_q)\right]=0$$ E8

Further, the normalized condition of phase volume void is written as in equation (8).

$${\displaystyle \sum _{q=1}^n{a}_q=1}$$ E9

Momentum equations:

VOF model takes all phases sharing same moving velocity, so the momentum equation is as follow.

$$\frac{\partial }{\partial t}(\rho \overrightarrow{\upsilon })+\nabla \cdot (\rho \overrightarrow{\upsilon }\overrightarrow{\upsilon })=-\nabla p+\nabla \cdot [\mu (\nabla \overrightarrow{\upsilon }+\nabla {\overrightarrow{\upsilon }}^T)]+\rho \overrightarrow{g}+\overrightarrow{F}$$ E10

In which, the density and viscosity are the weighted mean values calculated as follows:

$$\rho ={\displaystyle \sum a_q{\rho }_q}$$ E11

$$\mu ={\displaystyle \sum a_q{\mu }_q}$$ E12

F presents all the volumetric force except gravity, as described in equation (12). In which, σ is the tension force of water, and κ,α are the Surface curvature and liquid void at gas-liquid interface respectively., n presents outward unit normal vector of gas-liquid interface

$$\begin{array}{l}F={\alpha }_k(x)\alpha (x,t)\\ k=\frac{1}{\left|n\right|}[(\frac{n}{\left|n\right|}\cdot \nabla ]-\nabla \cdot n]\end{array}$$ E13

Turbulent modl equations:

Due to the rotation, the turbulent flow mostly is eddy flow, a modified turbulent model of RNG k-ε [10] in equation (13) and (14) are used in the work.

$$\frac{\partial }{\partial t}(\rho k)+\nabla \cdot (\rho k\overrightarrow{\upsilon })=\nabla \cdot ({\alpha }_k{\mu }_{eff}\nabla k)+G_k+G_b-\rho \epsilon -Y_M+S_k$$ E14

$$\frac{\partial }{\partial t}(\rho \epsilon )+\nabla \cdot (\rho \epsilon \overrightarrow{\upsilon })=\nabla \cdot ({\alpha }_{\epsilon }{\mu }_{eff}\nabla \epsilon )+c_{1\epsilon }\frac{\epsilon }{k}(G_k+c_{3\epsilon }{G}_b)-c_{2\epsilon }\frac{{\epsilon }^2}{k}-R_{\epsilon }+S_{\epsilon }$$ E15

$$c_{1\epsilon }=1.42,c_{2\epsilon }=1.68,{\alpha }_k={\alpha }_{\epsilon }\approx 1.39$$ E16

Where, G _k is turbulent kinetic energy by average velocity gradient, and Gb presents buoyancy production for turbulent kinetic energy. Y _M is contributed by the fluctuate expansion item for compressible flow, α _k and α _e are the reciprocal of effective Prandtl number for k and ε respectively. S _k and Sε are volumetric source items. Ρ and μ are weighted mean density and viscosity as defined previously in equations (10) and (11).

Incompressible flow with outer resource input is taken in this work, thus, G _b = 0, Y _M = 0 and S _k = S ε= 0.

The effective viscosity μ_eff in RNG k-ε model is related as in equation (13).

$$d\left(\frac{{\rho }^2}{\sqrt{\epsilon \mu }}\right)=1.72\frac{\stackrel{⌢}{\upsilon }}{\sqrt{{\stackrel{⌢}{\upsilon }}^3-1+C_{\upsilon }}}d\stackrel{⌢}{\upsilon }$$ E17

$$\stackrel{⌢}{\upsilon }=\frac{{\mu }_{eff}}{\mu },C_{\upsilon }\approx 10$$ E18

R_ε in equation (14) was given in equation (15).

$$R_{\epsilon }=\frac{C_{\mu }\rho {\eta }^3(1-\frac{\eta }{{\eta }_0}){\epsilon }^2}{1+\beta {\eta }^3}\frac{{\epsilon }^2}{k}$$ E19

$$\eta =S\frac{k}{\epsilon },{\eta }_0=4.38,\beta =0.01$$ E20

3. Geometry structure and calculation strategy

3.1. Geometry structure and mesh

The geometry structure of shaken flask is a unbaffled shake flask with a nominal volume of 250mL (Schott, Mainz, Germany), the same as Büchs used[9]. The shaker is orbital shakers with rotary diameter of 5cm. In order to decrease the calculation capacity, only the bottom half of flask with 5cm height as shown in figure 2, was considered because all the liquid is distributed in the zone and upper gas phase would not make difference for the liquid distribution when rotation.

Figure 2.

Schematic diagram of the simulation zone of the flask

The mesh of the calculation zone was plotted with Grid type by Gambit3.20. Grid independence study was preliminary carried out with three different mesh densities of close (49919), middle (29779) and coarse (12490), and found middle (29779) mesh cells was reasonable economic mesh good results. In order to enhance the convergence during iteration, smooth treatment of mesh was carried out to decrease the maximum skewness.

3.2 Dynamic mesh strategy

In CFD calculation, the quality of mesh would produce important effects on convergence and calculation efficiency. However, for those cases bearing moving boundary, the movement of the boundary would cause mesh cell stretched and skewed, worsen the quality of interior mesh cells related with moving area. In fluent, such problem was modified by DM technique. The principle of DM is based on local remeshing with size function strategy. In which, those cells of which the cell skewness or length scale is bigger than maximum cell skewness or maximum length scale, respectively would be labred and reunited first, and then remeshed to new required mesh cells according to size function factor γ. The size function factor γ is determined according to equation (17). More details of DM technique can be refereed in the help documents for Fluent 6.20[11].

$$\gamma =\{\begin{array}{l}1+m{d}_b^{1+2n},m0\\ 1+m{d}_b^{\frac{1}{1-n}},m0\end{array}$$ E21

3.3 Calculation conditions

The density and viscosity for both water and air are 1.225 kg m^-3, 1.7894×10^-5 Pa s, 1000 kg m^-3, 0.001 Pa s respectively. Surface tension of water is 0.072 N m.

In the simulation, water load of 15ml, 25ml and 35ml were investigated under different rotation speed of 100rpm, 150rpm, 200rpm and 250rpm respectively.

The equations were solved with 1st implicit with standard wall function and wall adhesion[12].The pressure interpolation scheme adopted was PRESTO (Pressure staggered option), which is useful for predicting highly swirling flow characteristics prevailing inside the flask. In order to reduce the effects of numerical diffusion, higher order discretization schemes are recommended. Accordingly, a third order accurate QUICK scheme was used for spatial discretization.

Fixed time step, Δt was used in the whole calculation, which was determined by equation (18).

$$\Delta t\le \frac{\Delta L}{u}$$ E22

Where, ΔL is the cell length of the cell connected moving boundary. It cab be estimated as 0.005m in this work. u is the maximum tangular velocity of shaken flask. As the results, $$\Delta t\le 0.003s$$ in this mesh system. Here $$\Delta t=0.001s$$ was adopted in all of the calculations.

The convergence criteria is 10^-3 in residual error. In our calculation, all the flow under different water load and rotation speed would arrive “steady state” in 2 second rotary time.

4. Results and discussion

4.1. Validation of the CFD results

The liquid shape and distribution in shaken flask with 25ml water load and rotation speed of both 250rpm (A) and 150rpm(B) by CFD simulation and by experimental and simplified model[9] calculated results were shown in figure 3. It can be observed that the simulated liquid distribution and shape have a good agreement with referred photographes[9] (A1,A2; B1,B2). The minor difference is possibly by the difference of viscosity in referred and this simulation. Bigger water viscosity could prevent the crawl of water during rotation, thus the maximum height of water crawl is a little smaller than this work.

4.2. The periodical position change of liquid phase

Figure 4 gives the periodical position change of liquid during rotation of 150rpm, 25ml water load. The same reference coordinate system was adopted for the nine angular positions. It can be found that as the flow arrived steady state, liquid in flask would periodically scan the relatively flask wall without any shape changes. It well agrees with the observation of flask movement.

4.3. Liquid phase shape

The rotation of shaken flask would produce centrifugal force to push water moving towards flask wall. Figure 5 gives the liquid shape under different rotation speeds with 25ml water load. Obviously, higher the rotation speed is, bigger the centrifugal force is and then smaller surface the water contact with the bottom surface of flask, and higher the water crawls to flask wall. Though the shaken flask is geometric symmetry, and the liquid phase position was changed periodically, the shape of water was not symmetrical. It is the combination effects of liquid rotary inertia and wall adhesion. Furthermore, with the increase of rotation speed, the gas-liquid interface would depressed deterioratedly, which means bigger area for gas-liquid contact.

Figure 3.

Simulated and referred liquid distribution in a 250mL shaken flask. Operating conditions: (A) V_L =25mL, N = 250 rpm, d₀ = 5 cm. (B) V_L =25mL, N = 150 rpm, d₀ = 5 cm. In which, A1, A2 and B1, B2 are from reference [8], A3 and B3 are the results of this work. It should be mentioned that the viscosity of referred data is 0.01 Pa s, ten times of the viscosity in this simulation.

Figure 4.

The liquid position changes in one rotation period of shaken flask with 25ml water load under 150rpm rotation speed.

Figure 5.

The liquid phase shape influenced by rotation speed with 25ml water load

4.4. Maximum height of liquid approached

As described previously, due to the rotation, liquid would be rejected outwards and crawled along flask wall. The effects of water load and rotation speeds on the maximum height of liquid approached are plotted in figure 6. It can be found that with the increase of both rotation speed and water load, the maximum height of liquid approached would also increase obviously. However, the increase of the maximum liquid height would dropoff due to the gravity limits when bigger rotation speed.

Figure 6.

Effect of liquid loading and rotation speed on the maximum height of liquid approached

4.5. Gas-liquid interfacial area

In shaken flask, the mass transfer is mainly dominated by the gas-liquid interfacial area. Bigger gas-liquid interfacial area always indicates good mass transfer, increasing rotation speed is a conventional method to improve mass transfer, such as oxygen supply for fermentation in shaken flask. Figure 7 presents the influences of both water load and rotation on gas-liquid interfacial area. It can be found, compared with water load, higher rotation speed does not lead to a bigger gas-liquid interfacial area, especially when the rotation speed is bigger than 150rpm. In fact, in most fermentation experiments in shaken flask, 150-200rpm rotation speed is a conventional choice. Too big rotation speed could not improve gas-liquid mass transfer obviously. More water load gives a rise of gas-liquid interfacial area, but it should be limited to the liquid crawls height on the wall in practice.

Figure 7.

Effects of rotation speed and water load on the gas-liquid interfacial area

4.6. Turbulent kinetic dissipation rate and volumetric power consumption

The rotation of liquid in shaken flask is produced by the power input from orbital shaking machine. How many the power consumed and how about the energy dissipated in flask are of interests to understand both mixing transfer and scale-up of bioprocess.

According to the turbulent model, turbulent intensity (I) is defined as in equation (19).

$$I=\frac{\upsilon \text{'}}{{\upsilon }_{avg}}$$ E23

In which, v’ is the.root mean square of turbulent fluctuation velocity, and v _avg presents the averaged flow velocity. The random fluctuation of micelles causes turbulence. So I indicates the degree of fluctuation and interaction between fluid micelles. Figure 8 gives the influences of rotation speed and water load on turbulent intensity. Higher rotation speed and more water load would almost linearly increase turbulent intensity, as the results, liquid mixing is better.

Figure 8.

Changes of turbulent intensity with different rotation speed and water load

Figure 9 plotted the changes of turbulent kinetic dissipation rates with different rotation speed and water load. It can be found that with the increase of rotation speed and water load, the averaged turbulent kinetic dissipation rates also increase greatly. Relatively, the influence of rotation speed is bigger on it. It is easy to understand that bigger rotation speed of shaken flask would input more power in it, such power mostly is dissipated by the collision between micelles in turbulent flow.

In most bioprocess, power consumption is usually taken as an effective criterion for process scaleup. Precisely understand of the power consumption in small apparatus or bench would help to improve the reliability of scaleup process and diminish pilot step. However, more details of power consumption of shaken flask are difficult, and in most cases, the power consumption in unbaffled Erlenmeyer flasks is calculated by empirical approach. Generally, the mechanical power introduced into the shake flask reactor during rotating motion is described as follows by Büchs et al. [13,14].

$$P=M(2\pi N)$$ E24

Where, P is the power input to flask, M is the net torque, and N is the shaken frequency (min^-1). Following the definition of the torque, the dimensionless power number for shake flasks is defined by Büchs et al. [14] as the modified Newton number:

$$Ne\text{'}=\frac{P}{\rho N^3{d}^4{V}_L^{1/3}}$$ E25

Figure 9.

Changes of turbulent kinetic dissipation rates with different rotation speed and liquid load

Moreover, the Reynolds number for shake flask is defined as

$$\mathrm{Re}=\frac{\rho N{d}^2}{\mu }$$ E26

and a relationship for dimensionless power in unbaffled shaken flasks is given by Büchs et al.[15]:

$$Ne\text{'}=70{\mathrm{Re}}^{-1}+25{\mathrm{Re}}^{-0.6}+1.5{\mathrm{Re}}^{-0.2}$$ E27

Thus, the volumetric power consumption in unbaffled shaken flask P/V_L (W/m³) can be calculated. On the other hand, according to the definition of turbulent kinetic dissipation rate ε, it means power consumption per unit mass. So $$\epsilon {\rho }_L$$ would also presents the volumetric power consumption rates The empirical approached results of P/V_L by equation (19-22) and averaged $$\epsilon {\rho }_L$$ by CFD simulation of all the cases of rotation speeds and water loads, is plotted in figure 11. It can be found that they both have a good linear relationship as in equation (20), and both the rotation speed and water load make no influence on the relationship.

$$\epsilon {\rho }_L=0.629\frac{P}{V_L}$$ E28

It means that the empirical value of volumetric power consumption P/V_L (W/m³) is much bigger than volumetric power consumption $$\epsilon {\rho }_L$$ (W/m³) by CFD. It should be noted that $$\epsilon {\rho }_L$$ is really dissipated by rotary liquid. Thus, $$\epsilon {\rho }_L$$ , or $$0.629\frac{P}{V_L}$$ by empirical approach is more practical as the power input criterion for scale-up of bioprocess directly from shaken flask.

Figure 10.

The relationship between volumetric power consumptions by empirical approach and CFD simulation.

5. Conclusion

A fluid dynamic model was prepared to calculate the flow behaviors in an Erlenmeyer shaken flask with a volume of 250mL moving in a shaker based on the combination of volume of fluid (VOF) model and RNG k-ε turbulence model. The numerical simulation was carried out using a commercial computational fluid dynamics software, Fluent, implementing the technique of Dynamic Mesh (DM). The effects of four different rotation speeds and three water load on the flow were investigated. The simulated results of liquid shape and distribution well agreed with experimental morphology in reference[9]. As the motion approached steady state, liquid would periodically traverse shaken flask wall with wall crawl driven by the centrifugal force and gravity without change of liquid shape and distribution. The shape of water was not symmetrical due the liquid rotary inertia and wall adhesion. The increase of both rotation speed and water load would give a rise of the maximum height of liquid approached Compared with water load, higher rotation speed does not lead to a bigger gas-liquid interfacial area, especially when the rotation speed is bigger than 150rpm. Bigger rotation speed and more water load would almost linearly increase the turbulent intensity to improve liquid mixing. With the increase of rotation speed and water load, the averaged turbulent kinetic dissipation rates also increase greatly. In order to under stand the power consumption in shaken flask, volumetric power consumption by empirical approach and CFD simulation were compared and it is found that there is good linear relationship between them without any influence by rotation speed and water load. And the empirical value of volumetric power consumption P/V_L (W/m³) is much bigger than volumetric power consumption $$\epsilon {\rho }_L$$ (W/m³) by CFD, probably indicating that $$\epsilon {\rho }_L$$ by CFD, or $$0.629\frac{P}{V_L}$$ by empirical approach is more practical as the power input criterion for scale-up of bioprocess directly from shaken flask.

A	–	Gas-liquid interface area, cm^2
a，b	–	Length in elliptical axis in x and y directions
H	–	Maximum liquid height approached
I	–	Turbulent intensity, %
k	–	Turbulent kinetic energy, m^2•s^-2
M	–	Net torque,
m	–	Mutation rate of size function
N	–	Rotation speed, rpm
Ne’	–	Modified Newton number
n	–	Mutation gradient of size function
P	–	Power input, W
Re	–	Reynolds number
t	–	Rotary time of shaken flask, s
Ux, Uy	–	Linear velocities of shaken flask in x and y directions, m•s^-1
u	–	Maximum linear velocities of shaken flask, m•s^-1
VL	–	Water load volume in flask, mL
v	–	Liquid flow velocity, m•s^-1
v’, vavg	–	Fluctuation and averaged velocity of liquid, m•s^-1
α	–	Volume void of fluid phases
γ	–	Size function factor
ε	–	Turbulent kinetic dissipation rate, m^2•s^-3
μ	–	Fluid viscosity of each phase weighted average, Pa•s
μeff	–	Effective viscosity, Pa•s
ρ	–	Fluid density of each phase weighted average, kg•m^-3
ω	–	Rotation angular velocity of flask, rad•s^-1
Subscripts
p	–	The p phase
q	–	The qphase

Table 1.