Instabilities due to turbulence through inlet jet in plunging jet bubble column

Author’s Accepted Manuscript Instabilities due to Turbulence THROUGH Inlet Jet IN Plunging Jet Bubble Column Ifsana Karim, Swapnil V. Ghatage, Mayur J. Sathe, Jyeshtharaj B. Joshi, Geoffrey M. Evans www.elsevier.com/locate/ces

PII: DOI: Reference:

S0009-2509(16)30260-3 http://dx.doi.org/10.1016/j.ces.2016.05.019 CES12958

To appear in: Chemical Engineering Science Received date: 23 November 2015 Revised date: 7 May 2016 Accepted date: 11 May 2016 Cite this article as: Ifsana Karim, Swapnil V. Ghatage, Mayur J. Sathe, Jyeshtharaj B. Joshi and Geoffrey M. Evans, Instabilities due to Turbulence THROUGH Inlet Jet IN Plunging Jet Bubble Column, Chemical Engineering Science, http://dx.doi.org/10.1016/j.ces.2016.05.019 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1

Instabilities due to Turbulence through Inlet Jet in Plunging Jet Bubble Column Ifsana Karim1, Swapnil V. Ghatage1, Mayur J. Sathe2, Jyeshtharaj B. Joshi3, Geoffrey M. Evans1* 1

Discipline of Chemical Engineering, University of Newcastle, Australia

2

Department of Chemical Engineering, Louisiana State University, USA

3

Homi Bhabha National Institute, Anushaktinagar, Mumbai 400 094, India

*

Correspondence to. Tel.: +61 240339068, Fax: +61 240339095, [email protected]

Abstract Bubble columns, where the liquid is in the continuous phase and the gas phase is on the form of dispersed bubbles, are widely utilised in many industrial applications. When designing these bubble columns it is important to be able to predict the conditions that mark the transition from homogeneous (bubbly) to heterogeneous (churn-turbulent) flow. The transition is a function of many variables, including: liquid and gas superficial velocities, bubble diameter, and liquid and gas physical properties. Linear stability analysis (LSA) has been successfully applied by many researchers to determine the gas volume fraction at which the instability takes place; and correspondingly a one dimensional stability factor, f1, has been proposed (see for example Joshi et al., 2001). Briefly, the LSA analysis utilises the velocity fluctuations of both the dispersed and continuous phases associated with the specific energy dissipation rate of the two phase mixture. Typically, the specific energy dissipation rate is correlated with the density of the bed. The previous analysis, however, does not consider the energy input which associated with the turbulence intensity (velocity fluctuations) of the incoming liquid stream. Usually, this component can be ignored in sparged bubble columns because its magnitude is relatively small and it also decays in the axial direction. In plunging liquid jet bubble columns the liquid is introduced as a high speed jet that entrains gas which is then broken into fine bubbles in the Mixing Zone. The bubby mixture then

2

passes into the Two Phase Flow Zone where instabilities can be generated. The Mixing Zone is a region of high energy dissipation resulting in relatively large liquid velocity fluctuations, which can directly influence the instability of the Two Phase Flow Zone. In this study the existing linear stability analysis is modified to include the influence of inlet liquid velocity fluctuations on the stability parameter, f1. The modified theory is applied to the previous work of Evans (1990) for a plunging liquid jet bubble column to determine the critical gas volume fraction at which transition takes place in the Two Phase Flow Zone. In order to apply the model, drift-flux analysis has been used to obtain bubble diameter as a function of gas and liquid superficial velocities, and computational fluid dynamics has been utilised to quantify the velocity fluctuations of the liquid exiting the Mixing Zone.

Keywords: Bubble column, regime transition, instability, turbulence, drift flux analysis.

Nomenclature A

parameter for stability criterion

(-)

B

parameter for stability criterion

(m/s)

C

parameter for stability criterion

(m2/s2)

CV0

virtual mass volume coefficient

(-)

D

column diameter

DG

gas phase dispersion coefficient

(m2/s)

DG

liquid phase dispersion coefficient

(m2/s)

DN

jet nozzle diameter

(m)

dB

bubble diameter

(m)

F

parameter for stability criterion

Fint.

interphase force exchange term (per unit volume)

(m)

(1/s) (kg/m2.s2)

3

f1

stability function defined by Eq. (2)

(-)

(f1)D

stability function denominator

(kg9/m20.s9)

(f1)N

stability function numerator

(kg9/m20.s9)

G

parameter for stability criterion

gz

gravity (has negative sign since acting in downward (negative) z direction)

k

turbulent kinetic energy

(m/s) (m/s2) (m2/s2)

K0,1,2 constant used in stability criterion

(-)

K3

parameter used in stability criterion, defined in Eq. (26)

(-)

Mo

Morton number (gµL4/( L23))

(-)

L

length

P

pressure

Re

Reynolds number (L VB∞dB/µL)

rB

bubble radius

t

time

(s)

Ta

Tadaki number (ReMo0.23)

(-)

UL

liquid superficial velocity

(m/s)

u

liquid interstitial velocity

(m/s)

u’

liquid fluctuating velocity in the liquid surrounding the bubbles

(m/s)

(m) (kg/m.s2) (-) (m)

u’inlet liquid fluctuating velocity of the inlet liquid

(m/s)

VG

gas superficial velocity

(m/s)

VG’

gas drift-flux velocity

(m/s)

Vj

jet velocity exiting the nozzle

(m/s)

VB∞

gas (bubble) terminal velocity in an infinite liquid

(m/s)

VS

slip velocity (v-u)

(m/s)

v

gas interstitial velocity

(m/s)

4

v’

gas fluctuating velocity

Z

parameter of the stability criterion

z

axial (vertical) position

(m/s) (m2/s2) (m)

Greek Letters α, β

fitted parameters in Eq. (35)

(-)

Ф1, Ф2 fitted parameters in Eq. (37)

(-)

ϵG

gas (dispersed phase) volume fraction

(-)

ε

specific energy dissipation rate

(m2/s3)

G

gas (dispersed phase) density

(kg/m3)

L

liquid (continuous phase) density

(kg/m3)

µL

liquid (continuous phase) dynamic viscosity

(kg/ms)



liquid surface tension

(kg/s2)

Subscripts G

gas (dispersed) phase

L

liquid (continuous) phase

inlet

inlet to Two Phase Flow Zone

L

liquid (continuous) phase

D

denominator

MZ

Mixing Zone

N

numerator

1.

Introduction Efficient phase interactions are central to all multiphase reactors, especially for gas-liquid

operations carried out in bubble columns, stirred tanks, etc. Bubble columns are often favoured over other gas-liquid multiphase equipment due to their high efficiency and simplicity of construction.

5

Typically, bubble columns have an upflowing gas phase whilst the liquid phase is either stationary or downflowing. Numerous variations exist, however, depending on the operation being performed. For example, plunging liquid jet bubble columns are used for mineral flotation processes as they produce very fine bubbles and have high gas volume fractions—where both of these characteristics are beneficial for the recovery of the valuable mineral product. In the plunging jet bubble column the gas and liquid phases enter through the top of the column. The liquid phase is injected through the centre in the form of a jet, which entrains the gas as bubbles. The bubbles travel in a downward direction against gravity due to the drag force exerted by the liquid. Two zones with different flow patterns are generated within the column. The Mixing Zone is at the top of the column where the liquid jet plunges into the fluid below and creates a highly circulating two-phase flow with a high rate of energy dissipation. The Mixing Zone performs two important tasks, namely: (a) entrainment of gas into the liquid, and (b) breakup of larger bubbles into finer bubbles, thus generating the gas-liquid dispersion. The dispersion exits the Mixing Zone and into the Pipe Flow Zone before being discharged from the outlet located at the bottom of the column. As the name implies the Pipe Flow Zone has the typical characteristics of two phase flow inside a pipe, i.e. uniform pressure drop in the axial direction. The exchange of momentum between the gas and liquid phases in plunging jet column is a complex phenomenon. Bubble columns can operate in either homogeneous or heterogeneous regime. In the homogeneous regime, the gas volume fraction and liquid velocity are deemed to be uniform over the cross-section of the column. Typically, either no or very little bubble breakup or coalescence takes place. In the heterogeneous regime, both the gas volume fraction and velocity profile of the liquid phase become non-uniform and the bubble size usually increases significantly due to coalescence. Heterogeneous flow is associated with higher turbulence and generally provides better heat, mass and momentum transfer rates. Conversely, it can be detrimental to the efficiency of some processes such as flotation due to the disruption of bubble surfaces (interfacial area) onto which the recovered minerals have become attached to. The preference for either homogeneous or

6

heterogeneous regimes will depend on the particular application, and in order to optimise the system it is important to know the operating ranges for each of these regimes and the point at which the transition takes place. For a plunging liquid jet bubble column with a given geometry and liquid jet velocity regime transition is tied closely to the gas superficial velocity. At very low gas superficial velocities the operation regime is almost always homogeneous, whilst at much higher gas superficial velocities the operation regime will be heterogeneous. The homogeneous to heterogeneous regime transition occurs at some gas superficial velocity between these two limits. The regime transition for the plunging liquid jet bubble column is determined by the balance between stabilising and disturbing forces. In the homogeneous regime the stabilising forces dominate, suppressing the instability caused by perturbations such as fluctuation of velocity and volume fraction of either phase. In the heterogeneous regime the stabilising forces are weak, causing flow perturbations to grow. Such phenomena have been investigated using linear stability analysis. Briefly, perturbations are applied to the constitutive equations for two-phase flow, and the growth or decay of such perturbations are used as the indicators of the regime of operation (Joshi et al., 2001). The approach has led to the development of a stability function, f1, whereby if f1 is positive then the system is homogeneous; whilst if f1 is negative then the system is heterogeneous. The root(s) of this function (f1 = 0) with respect to the perturbation variable, which is usually the dispersed phase volume fraction, therefore represent the critical condition at which transition takes place. It has been observed that the stability function, f1, can exhibit multiple transitions (roots) from the homogeneous-to-heterogeneous regime and vice-versa. For example, Gibilaro et al. (1986) found that for a solid-liquid fluidized bed the system went initially (at low solids volume fraction) from a homogeneous-to-heterogeneous regime, but then underwent a second regime transition (heterogeneous-to-homogeneous) at a much higher volume fraction of around 0.8-0.85. Other examples of multiple root behavior for fluidised beds are given by Sundaresan (2003). Ghatage et al. (2014) found a similar behavior for bubble columns, but the difference was that the transition from heterogeneous-back to-homogeneous regime occurred at a critical gas volume fraction that was very

7

much closer (with a few percent) to the original homogeneous-to-heterogeneous transition. The suggestion that was proposed for this behavior was that the linear stability approach no longer applied once the first transition had occurred, and that the reverse transition point should therefore be ignored. The breakdown in the stability approach was due to both changes in bubble size due to coalescence/breakup, and the flow no longer being one-dimensional as is assumed in the theory. In this study linear stability analysis has been applied to a plunging liquid jet bubble column, with the principal goal of being able to predict when regime transition takes place in the Pipe Flow Zone as a function of gas input for a given column geometry, fluid properties, and jet velocity. In addition, the study aims to address the following points that remain uncertain. With respect to the plunging liquid jet bubble column currently being considered, these factors include: 1.

Liquid phase velocity fluctuations generated by the energy dissipation of the plunging liquid jet within the Mixing Zone.

2.

Change in bubble size due to coalescence/breakup and its meaning in terms of stability.

3.

Significance of each of the roots that constitute the stability criteria.

In order to explore each of these factors the experimental data of Evans (1990) is utilised, coupled with drift-flux analysis and computational fluid dynamics to obtain estimates of bubble size and velocity fluctuation of the liquid phase entering the Pipe Flow Zone, respectively.

2.

Modelling Approach

1.1

Linear stability parameter numerator and denominator analysis Jackson (1963) first applied linear stability analysis to multiphase flows, whereby small

perturbations are introduced into the governing equations of continuity and motion for each of the phases. Since then there have been numerous studies on the topic (e.g. Anderson and Jackson, 1969; Lahey, 1991; Shaikh and Al-Dahhan, 2007; Nedeltchev and Shaikh, 2013; Ghatage et al.; 2014). Of interest to this study is the work of Joshi et al. (2001) who developed a stability criterion for both

8

upflow and downflow bubble columns in either batch or continuous operating modes. They proposed the following stability criteria marking the transition from homogeneous to heterogeneous flow:

f1  0 ,

(1)

where

 A  G F   B 2 A  Z  C   B2 4  A  G F   B 2  f1  1   A  Z  C   B2 4 A  Z  C   B2 4 2

2

.

(2)

In Eq. (2), for a gas-liquid bubbly flow the parameters A, B, C, F, G, Z are functions of ϵG, ρL, ρG, gZ, dB, VB∞, VG, UL, CV0, K0, K1, K2, K3. Expressions for these terms are reported in the literature (e.g. Joshi et al., 2001). As mentioned previously, a system is deemed unstable when f1 = 0; and a typical plot for bubble column of f1 versus ϵG is shown in Fig. 1. Many studies based on linear stability approach showed presence of two transition points, however, only stable-unstable transition can be seen practically. It can be seen from Fig. 1 that f1 crosses the horizontal axis at point 1 and 2; where it is supposed that the point 1 marks the transition from stable to unstable condition as ϵG is increased. Similarly, for point 2 the interpretation is that the system transition from unstable to stable condition as ϵG is further increased. Often, the separation between points 1 and 2 is very small, and the physical meaning of the transition from stable-unstable-stable with little increase in ϵG is not immediately apparent. The argument often used to that once the system becomes unstable (point 1) the model equations are no longer valid and therefore the (mathematical) curve to the right point 1 no longer applies in real case. This seems to be quite a reasonable argument. In this study a different approach to the analysis and interpretation of the stability criteria is proposed. Firstly, it is assumed that the expressions for the parameters A, B, C, F, G, and Z can be expressed in terms of their denominators (denoted by D) and numerators (denoted by N), e.g. A=AN/AD. Using the relevant expression given in the literature then:

9

G K0 1 G    L K0 1 G    L CV 0 1  2 G    L K0 1 G  A ; A N  2 AD L K0 1 G 

(3)

 2G K 0VG 1 G 3  2LVG CV 0 1  2 G 1 G 2  2L K 0U L G2 1 G   BN  B   3 BD    L K 0 G 1 G  ;  2 L CV 0U L G2 1  2 G     3   L K 0 G 1 G  

(4)

 G K 0VG2 1 G 4   L CV 0VG2 1  2 G 1 G 3   L K 0U L2 G3 1 G   CN  C   4 CD    L K 0 G2 1 G  ;   L CV 0U L2 G3 1  2 G     4   L K 0 G2 1 G  

(5)

2

F

2

g Z G  G  L  FN  ; FD L VG 1 G   U L G 

(6)

2 2   GN g Z  G  L  VG 1 G   U L G  G  ; GD L 1 G  VG 1 G   U L G 

(7)

 K  1  3  0 G G     g K d V 1    U    Z Z 2 B G G L G 3  Z  N    C 1  K  1  2   .     V 0 3 G G 3 ZD  K 0 1 G  VG 1 G   U L G     2 2    K K  1    0 3 G  G

(8)

Substituting Eqs. (3)–(8) into Eq. (2) gives:

 f1  N f1   f1 D



Z D CD  2 AN BD GN FD  AD BN GD FN 

2

AD GD2 FN2  4 AN BD2  CD Z N  CN Z D   AD BN2 CD Z D 

.

(9)

Note that the expressions where K0, K1, and K2 appear are given by Bhole and Joshi (2005) as:

 3  2 1 G   CV  CV 0   K 1    , G  0 

(10)

CV 0  0.5  K1 Ta  ,

(11)

DG  K2 d B G VS ,

(12)

10

where DG is the gas phase dispersion coefficient, VS is the slip velocity, Ta is the Tadaki number. Generally reported values for K0, K1, and K2 are 1, 0.143 and 3, respectively; and are the values adopted for this study. It can be seen from Eq. (9) that the stability parameter, f1, can be expressed in terms of its numerator and denominator functions (f1)N and (f1)D, respectively. The reason for reporting f1 this way is that it can be used to highlight particular features of the stability parameter curve as it changes as a function of gas volume fraction. To assist with the highlight, both (f1)N and (f1)D are each expressed as polynomial functions of the gas volume fraction. Using the previous work of Ghatage et al. (2014), and following extensive manipulation of the various terms, (f1)N and (f1)D can be expressed as following polynomials:

( f1 ) N  Y1 Y2Y3  Y4  2 ,

(13)

where: Y1  23L  G  L  g z G3 1 G  VG 1 G   U L G  K0 ;

(14)

2 Y2  VG 1 G   U L G2  ;  

(15)

2 Y3  G 1 G   L 1 G  G L CV 0 1  2 G  ;  

(16)

GVG 1 G 3  L CV 0VG 1  2 G 1 G 2 Y4   LU L G2 1 G   L CV 0U L G2 1  2 G 

(17)

2

10

 : 

Similarly,

( f1 ) D  X1 4 AN  CD X 2 X 3  CN X 4   X 5  ,

(18)

where: X1  45L  G  L  g z2 G4 1 G  VG 1 G   U L G  ;

(19)

X 2   g z K2 d B VG 1 G   U L G ;

(20)

3 2 X 3 G  K0 1 G   CV 0 1  K3  G2 1  2 G   K3 K0 G 1 G   ;  

(21)

2

13

3

11

X 4  K0 1 G  VG 1 G   U L G  ;

(22)

X 5  K0 1 G  VG 1 G   U L G  .

(23)

3

3

Equations (13)-(23) show that both (f1)N and (f1)D can be expressed as polynomial functions that, individually, are continuous have no discontinuities. It can be seen that the terms Y1 and X1 that appear in (f1)N and (f1)D, respectively, will result in multiple roots at ϵG = 0 and ϵG = 1 in both the numerator and denominator. Hence, the resultant f1 curve will have a zero (root) value at both of these gas volume fractions, which is perhaps not surprising given that volume fraction outside of the range (0 ≤ϵG ≤1) don’t exist in real systems. It is also interesting to note that both the Y1 and X1 terms will result in an additional root when: 0  VG 1 G   U L G  ; 3

(24)

that occurs at the no-slip gas condition, i.e.:

G noslip 

VG . VG  U L

(25)

When either (f1)N or (f1)D is plotted as a function of ϵG this root just touches the ϵG axis, which is different to what occurs at the roots at ϵG equal to 0 and 1, where f1 actually crosses the horizontal axis. The behavior of the real system for both types of roots is discussed later in this paper. Equations (13)-(23) show that both (f1)N and (f1)D, include an un-factorised component, given by [Y2Y3-Y4]2 and [4AN(CDX2X3-CNX4)+X5], respectively. When the relevant substitutions are made, which include Eqs. (15)-(17), (3), (5), (20)-(23), the resultant equations are also polynomials of order at least 12. Unlike the cases for the roots previously identified at ϵG equal to 0, 1 and no slip condition, VG/( VG + UL), the root(s) that might be generated from the polynomial expressions for [Y2Y3-Y4]2 and [4AN(CDX2X3-CNX4)+X5] are functions of the system operating conditions that include VB, VB, L, G, , gz, L and K0-3. Arguably, therefore, it is these roots that will provide greater

12

insight into the stability of the system as a function of conditions that go beyond just gas and liquid superficial velocity and where slip between the phases is present. 1.2

Bubble diameter determination As indicated above the critical gas volume fraction—apart from at the no slip condition—that

marks the onset of instability between the realistic operating range of 0<ϵG<1 can be obtained by finding the roots of [Y2Y3-Y4]2 and [4AN(CDX2X3-CNX4)+X5]. To do so, requires values for ρL, ρG, µL, , gz, dB, VG and UL; and for a given system these values are known. Also, CV0 is given by Eq. (11), an appropriate expression for VB∞ can be obtained from the literature (e.g. Clift et al., 1978), and numerical values for K0, K1 and K2 are reported above. The remaining parameters that need to be determined, therefore, are dB and K3. In many gas-liquid bubble column operations the input bubble diameter, dB, is a variable of the system depending on the gas injection method. For example, for orifices, porous spargers, venture devices, impinging gas jets, plunging liquid jets, etc., the resultant bubble size will almost certainly be a function of the gas flow rate. It is only under very controlled conditions that the bubble size can be kept constant whilst the gas flow rate is varied. It is for these reasons that the bubble size is usually measured experimentally, but the exercise becomes very difficult at gas volume fractions above about 0.2—which is the region where the instability is most likely to take place. It is possible, however, to obtain a numerical value for the bubble diameter by applying one-dimensional drift-flux analysis to the operating conditions and the measured gas volume fraction. This approach was used by Evans (1990) for a plunging liquid jet bubble column, and is re-applied in this study to obtain the corresponding bubble diameter estimations needed for the stability analysis. 1.3

K3 determination The only remaining variable that needs to be quantified before the stability analysis is the

liquid velocity fluctuation, u’, which unlike the fluctuating velocity of the gas (bubble) phase, v’, has no generally accepted expression reported in the literature. Typically, u’ is taken to be equal to zero,

13

which is perhaps a reasonable assumption for bubble columns operated (1) in batch mode, where there is no net liquid flow; or (2) with a gas sparger where the velocity fluctuation (turbulence intensity) of the inlet liquid phase, u’inlet, might be very low. For plunging liquid jet bubble columns, the high speed liquid jet generates significant liquid velocity fluctuations at the top of the column in the Mixing Zone. These liquid velocity fluctuations are maintained at the lower boundary of the Mixing Zone due to the constant action of the plunging liquid jet. It is therefore not appropriate to assume that u’inlet entering the Two Phase Flow Zone below is equal to zero. The mean liquid velocity fluctuation in the Two Phase Flow Zone, v’, is the combination of the liquid fluctuations that occur within: (1) the virtual mass volume, CV0, of liquid that is influenced by the fluctuating gas phase; and (2) the bulk liquid—defined as the volume of liquid not within the virtual mass volume—that is influenced only by the inlet liquid velocity fluctuation. The corresponding liquid fluctuating velocity profile is illustrated in Fig. 2, from which an average value can be obtained by integrating over the entire (bulk plus virtual mass) liquid volume. It’s worth noting at this point that the liquid velocity fluctuation in the virtual mass volume will remain constant along the Two Phase Flow Zone since the fluctuation of the bubbles is continuously providing energy input into the liquid. The liquid velocity fluctuation in the bulk liquid will dissipate along the axial direction, however, since the source of their existence, namely the liquid jet plunging into the Mixing Zone, is no longer being applied. In a previous study (Ghatage et al., 2014), it was proposed that the u’inlet exiting from the Mixing Zone of a plunging liquid jet bubble column can be expressed as some multiple, K3, of the fluctuating velocity of the gas (bubble) phase, v’, within the Two Phase Flow Zone, i.e.:

u 'inlet  K3v ' .

(26)

The numerical value for K3 in Eq. (26) will be determined by the turbulence generating mechanisms/devices being applied at the inlet; and for the plunging liquid jet, K3 will be a function of column and jet diameter, jet velocity, gas flow rate and liquid and gas physical properties.

14

Many computational fluid dynamics (CFD) studies reported the hydrodynamic behaviour of bubble columns, but relatively few have actually focused on the mixing and transition behaviour taking place between the gas and liquid phases (e.g. Joshi et al., 2001; Ekambara et al., 2005; Monahan et al., 2005; Shaikh and Al-Dahhan, 2007; Upadhyay et al., 2009). Ekambara et al. (2005) carried out CFD simulation of bubble columns in 1D, 2D and 3D geometry and studied mixing behavior and axial dispersion coefficient using k-ε model. The authors found that their 1D and 2D models over-predicted by 20-100 and 50-200 percent, respectively; with the 3D model being most promising especially in terms of eddy diffusivity. Upadhyay et al. (2009) investigated the hydrodynamics of a downflow bubble column and the effect of phase throughputs on profiles of gas volume fraction and intermixing behaviors. Their study involved experimental as well as numerical analysis. They performed Eulerian-Eulerian simulations in a 2D axisymmetric geometry with standard k-ε model to visualize the turbulence properties in the column. They considered drag, lift and virtual mass force and recommended that the CFD model was acceptable to predict average cross-sectional volume fraction values, however, unable to predict the phase distribution accurately. In this study, the computed CFD simulation outputs of liquid phase turbulent kinetic energy and gas phase volume fraction been used to determine the relationship between K3 and the system conditions in order to apply the stability analysis methodology.

3.

Model Validation

3.1.

Experimental Measurements The modelling analysis has been applied to the experimental work of Evans (1990) for a

plunging liquid jet bubble column. Briefly, the experimental system consisted of a vertical column (0.044 m diameter, 1.000 m long) which housed a plunging liquid jet (0.0047 m diameter, 11.5 m/s velocity). Gas (air) was introduced in the enclosed headspace and was entrained by the liquid (water) jet to create a downflowing bubble mixture. Pressure tappings at the base of the column, at 0.950 and

15

1.000 m, were used to measure the gas volume fraction, ϵG, as a function of column liquid, VL, and gas, VG, superficial velocities. The superficial liquid velocity was fixed at 0.134 m/s, whilst the gas velocities ranged from 0.028-0.137 m/s. The experimental data for ϵG (closed circles) versus VG is shown in Fig. 3. Also shown in the figure is the gas volume fraction based on a no slip condition, (ϵG)noslip, and computed from Eq. (25). As expected, it can be seen that ϵG increased with increasing VG and that the measured ϵG data points fell above the (ϵG)noslip curve as there was slip between the gas and liquid phases due to the buoyancy of the bubbles. Relevant to this study are two points. Firstly, the curve for the ϵG measurements appears to follow a continuous trend, such as that of a power law, with no obvious deviations that would suggest a transition from homogeneous (stable) to heterogeneous (unstable) conditions. Secondly, according to Eq. (14) the system will become unstable, as defined by (f1)N being equal to zero, once the measured ϵG becomes higher than the corresponding (ϵG)noslip value. According to the information presented in Fig. 3 this condition was reached for all VG values, and consequently the flow should have been heterogeneous (unstable). Visual observation clearly indicated that for the smaller superficial gas velocities this was not the case since the gas phase was homogeneously distributed as very fine bubbles within the liquid phase. Therefore, it seems that it is not an absolute certainty the system will become unstable once the measured ϵG goes above the no slip condition—at least not for the cocurrent, gas-liquid, downflow case. The reason for such an apparent violation of the mathematical model can be found by consideration of how the bubble size changes as a function of the conditions inside the column. 3.2.

Bubble size estimation An estimate of the bubble diameter, dB, is required in order to apply the stability analysis.

The bubble diameter has been obtained by application of drift-flux analysis to the measured ϵG values for each VG setting. An application of the methodology is shown in Fig. 4. Firstly, the drift-

16

flux velocity, VG’, for each data point (closed circles) has been calculated according to the following expression for the drift-flux operating line:

V 'G  VG 1 G  U L G .

(27)

Secondly, the drift-flux curve (solid line) has been calculated (see Evans, 1990), where the bubble radius has been varied until the drift-flux curve and operating line intersect with each other at the measured gas volume fraction. For the example given, the conditions at the intersection (open circle) are: VG = -0.031 m/s, UL = -0.134 m/s, and (ϵG)measured = 0.20; resulting in a calculated bubble radius of rB = 0.000157 m. By plotting different drift flux curves for each of the experimental observations, rB values have been obtained as a function of gas superficial velocity, and these are plotted in Fig. 5 along with the corresponding gas volume fraction, ϵG. It can be seen that rB increases with increased VG, which can be accounted for by consideration of the energy dissipation rate as well as bubble coalescence and breakup events taking place within the Mixing Zone at the top of the plunging jet bubble column. The values of rB are less than 0.000800 m even at the highest reported value of VG. Such fine bubbles are usually present only in homogeneous two phase flow, and subsequently it’s not apparent at what gas input, if at all, the system has become unstable. The bubble radius estimates do indicate, however, a significant increase in size at a VG of around 0.07 m/s; and possibly this change in behavior is the result of an instability being reached within the system. Before this possibility can be explored by applying the stability analysis, a quantitative value of the velocity fluctuation of the incoming liquid phase, u’inlet, is required. To do this, in accordance with Eq. (26) an estimate of K3 is required. 3.3.

K3 estimation

3.3.1 Computational fluid dynamics modeling An estimate of K3 for the experimental operating conditions was obtained by application of the computational fluid dynamics package, FLUENT-14. Briefly, an Eulerian-Eulerian modelling approach was applied given that the gas phase volume fraction was expected to vary widely

17

throughout the flow domain. Standard continuity and momentum equations were applied to both the gas and liquid phases. For example, the relevant equations for the liquid and gas phases in three dimensional geometry are:

L

G

L

 1 G  t





 L 1 G  uz  0 ,

(28)

 G  G  G vz  0 , t





 1 G  u z t



(29)

    1   u u    L

G

z

z

 1 G   PL

   

  1 G   L u z  u z 

T

 , 

(30)

 1 G   L g z  Fint.

G

 G vz  G  G vz vz  t





   

 G  PG   G G vz  vz 

T

  ,

(31)

 G G g z  Fint. where P is pressure, µ is absolute viscosity and Fint. denotes the interphase exchange force. The standard k-ε turbulence model (Launder and Spalding, 1974) was used for the simulations as it is accurate and robust. The drag force between the gas and liquid phases was modelled using the Schiller-Naumann drag law. As indicated previously, the computed bubble size was very small (dB<500 µm), so acceleration and virtual mass force effects were assumed to be negligible. All other forces, e.g. Bassett, have also been neglected; and additionally the kinetic theory of granular flow was not applied given that the bubbles can coalescence and/or breakup. The computational geometry used for the plunging liquid jet bubble column is shown in Fig. 6. The column itself was 0.044 m diameter and 1.000 m long. The liquid inlet at the top of the column was a 4.76 mm diameter nozzle (jet) facing downwards. The gas inlet was an annular opening at top bounded by the outside diameter of a nozzle and the inside diameter of the column. The entire volume of the column was discretized using 3D structured mesh and a total of around 0.3

18

million nodes were used. As shown in Fig. 6 the mesh coarsened outwards in the radial (x) direction, and along the axial (z) direction from the (top) inlet to (bottom) outlet. A uniform liquid velocity of 11.5 m/s was specified at the liquid inlet (nozzle), corresponding to liquid superficial velocity, VL, of 0.134 m/s inside the column. The turbulent intensity at the jet inlet was defined to be 5 per cent. In addition, the hydraulic diameter was specified as 0.0048 m. Constant pressure outlet (101,325 Pa) was defined at the bottom outlet and a no slip boundary condition was employed at all the walls. Simulations were carried out at three different superficial velocities of gas, VG, of 0.028, 0.049 and 0.069 m/s. The computation was undertaken for a full 3D unsteady system. All terms of the governing equations for unsteady state were discretised using the second-order upwind differencing scheme. In all the simulations, a time step of 1×10−4 s was chosen. The SIMPLE algorithm was employed for the pressure-velocity coupling. The convergence criterion (sum of normalized residuals) of 1×10−4 was set for all the equations. Initially the column was filled with liquid. At time t=0 s, liquid jet and gas were introduced in the column. The flow profiles took around 8-10 s to develop. The simulations were performed for 20 seconds and averaged data over 3 seconds are reported in the present study. Steady state was defined when there was less than 5 percent change in phase volume fractions. The computed outputs were gas phase volume fraction, liquid and gas velocities, pressure, turbulent kinetic energy, k, and specific energy dissipation rate, ε. The contour plots of gas volume fraction at gas superficial velocities of (A) 0.028, (B) 0.049 and (C) 0.069 m/s are shown in Fig. 7. It can be seen that the gas occupies the whole section above the liquid inlet. This portion of the column is known as the headspace where the liquid jet plunges before entering the bubbly mixture below. The gas volume fraction at the centre of the column is less than 10 per cent mainly because the high speed of the jet is providing less time for the gas to mix at the centre. As the gas-liquid mixture moves downward an annular recirculation region is observed between the expanding submerged jet and the column wall. The two phase mixture in this region has the highest volume fraction of the gas phase. This is mainly because gas gets entrapped in the

19

section, mixing slowly with liquid flowing downstream. The recirculation region below the free surface is known as the Mixing Zone, and is marked by strong recirculation, high energy dissipation rate, and is where the bubble size is largely determined that enters into the relatively quiescent Two Phase Flow Zone below. The gas volume fraction contour plots shown in Fig. 7 have been averaged over time and space across the flow area in the axial direction to obtain the mean gas volume fraction as a function of axial distance. The results are shown in Fig. 8 for the three gas superficial velocities, VG, of 0.028, 0.049 and 0.069 m/s. It can be seen that as expected G increased with increasing VG. For all three velocities there is a region of elevated gas volume fraction at the top of the column that marks the boundary of the Mixing Zone where recirculation occurs and the motion of the gas and liquid phases are directly influenced by the plunging liquid jet. Beyond the Mixing Zone is the Two Phase Flow Zone, where the simulations show that gas volume fraction steadily decreases in the axial direction. The measured gas volume fractions, monitored between 0.950-1.000 m, are also shown in Fig. 8 (large open circles); and it can be seen that the CFD predictions are in good agreement with the measurements. 3.3.2 Identification of inlet location As mentioned previously, the CFD simulations have been undertaken in order to obtain an estimate of the K3 term so that a stability analysis can be carried out. The positive comparison between the simulated and measured gas volume fractions gives confidence in proceeding with the utilisation of the simulated output values for determination of K3. As defined in Eq. (26), K3 is simply the ratio of u’inlet divided by v’, where both of these terms can be computed from the simulation output. For example, u’inlet, which is the liquid inlet fluctuating velocity, is equal to the square root of turbulent kinetic energy, k, at the liquid inlet, i.e.:

  kinlet , uinlet

(32)

20

where the “inlet” position is the end of the Mixing Zone since the stability analysis is being carried out for the flow inside the Two Phase Flow Zone. Similarly, the fluctuating gas phase velocity, v’, within the Two Phase Flow Zone can also be obtained from the simulated results and assumption that (Joshi et al., 2001):

v ' G U L  VG  .

(33)

Substituting Eqs. (32)-(33) into Eq. (26) and rearranging gives the following expression for K3 to be used for the Two Phase Flow Zone:

  kinlet K3   .  G U L  VG   two phase zone

(34)

Before proceeding with plotting K3 as a function of axial position, it is timely to present the axial profile for the turbulent kinetic energy, k, which is shown in Fig. 9. It can be seen that for all three gas superficial velocities, that towards the top of the column k increases rapidly to a peak value, corresponding to the centre of the recirculation region in the Mixing Zone, before rapidly decaying to some small baseline value—approximately 0.001, 0.002 and 0.004 m2/s2 for VG equal to 0.028, 0.049 and 0.069 m/s, respectively. The shape of the k curve is important for the K3 analysis, as it is the value at the inlet to the Two Phase Flow Zone, kinlet, that is required according to Eq. (34). From the graph shown in Fig. 9 it would appear that the axial positions for which the kinlet value should be computed are approximately 0.3, 0.6 and 0.7 m for VG equal to 0.028, 0.049 and 0.069 m/s, respectively. Further validation that it is reasonable to use these axial positions for the determination of K3, is given in Fig. 10, where the computed specific energy dissipation rate, , axial profile is presented. It can be seen that similar to the curves for the kinetic energy,  is highest in the Mixing Zone but decays rapidly upon entering the Two Phase Flow Zone as indicated by the heavy black lines given by the following relationship:





 z  LMZ 



,

(35)

21

where LMZ is the axial length of the Mixing Zone (0.3, 0.6, 0.7 m for VG equal to 0.028, 0.049 and 0.069 m/s, respectively) and α and β are fitted parameters. The numerical value for β was 1.0 for all three gas superficial values, whilst values for α were 0.013, 0.014 and 0.023 for VG equal to 0.028, 0.049 and 0.069 m/s, respectively. 3.3.4 K3 versus superficial gas velocity The axial profile for K3 is shown in Fig. 11, where K3 has been computed using the local averaged values for k and G at that axial location. It can be seen that, similar to the curve for the specific energy dissipation rate, the K3 curve reaches a peak value and then gradually decays along the axis of the column. Whilst it is reasonably clear that the boundary between the Mixing and Two Phase Flow Zones lies at 0.3 m for the lowest superficial gas velocity, it is not as obvious where the boundary is for the two higher velocities. The boundaries are much more easily identified from plots of k and  given in Figs. 9 and 10, respectively. The numerical values for K3 at the Mixing Zone exit determined from the graphs are 3.6, 2.9 and 4.9 for VG equal to 0.028, 0.049 and 0.069 m/s, respectively. The result that K3 varies with VG is not surprising given that it is dependent on the fluctuating gas velocity and will almost certainly be influenced by the energy dissipation in the Mixing Zone—which is itself a function of the gas input. The CFD simulations, in theory at least, will have accounted for all of the process influences on K3. However, the CFD simulations require additional processing to give numerical values for K3 that can be used in conjunction with the linear stability analysis. For this reason, an attempt has been made to correlate the three (simulated) K3 values with superficial gas velocity, VG, which was the process parameter varied as part of this study. The correlation required consideration of the boundary conditions under which K3 might be related to VG. Firstly, at no gas input into the system, i.e. VG = 0 m/s, then v’ will also be zero. The incoming liquid, however, can have some velocity fluctuation, i.e. u’inlet ≠ 0 m/s. From Eq. (26) it follows that at v’ = 0, then K3 → ∞. A functional form involving VG that will satisfy this relationship is:

22

K3 

1 . VG

(36)

where 1 is a constant. Equation (36) satisfies the condition at zero gas superficial velocity, where K3 is at (or at least approaches) infinity, and then decays as VG is increased until it eventually approaches zero at infinite VG. However, in a real bubbly/heterogeneous flow system the bubbles will become so densely packed that their fluctuations are reduced to zero, well before VG approaches infinity. At this superficial gas velocity limit, (VG)limit, then v’ will become zero, However, the same condition is not applied to the incoming liquid since it has not entered the densely packed structure. Therefore, at (VG)limit, then u’ ≠ 0 m/s and K3 → ∞. Hence, an additional term is required in Eq. (36) to satisfy this condition; and following a similar approach to that used for VG = 0 m/s, the resultant expression is:

K3 

1 2  . VG VG limit  VG

(37)

Equation (37) has been applied to the computed K3 values as shown in Fig. 12, where the fitted values for Ф1, Ф2 and (VG)limit are 0.079, 0.039 and 0.079 m/s, respectively. Whilst it is not relevant to comment on the agreement of the fit since three fitting parameters have been used to fit three data points, it is worthwhile noting that the trend of the minima for the three data points is able to be modelled using the functional form given by Eq. (37) using realistic values of the fitting parameters, i.e. (VG)limit of 0.079 m/s, which corresponded to visual observations that the column had become unstable at this velocity. Equation (37) has been applied using the fitted parameters to the three simulated data points to the full range of experimental superficial gas velocities (0.028≤ VG≤0.137 m/s) reported in Fig. 3. The resultant graph is shown in Fig. 13, where a nominal cut-off value of 50 for K3 has been chosen for VG beyond (VG)limit. Note that any value above about 20 will yield close to the same VG value for the onset of instability as discussed below. 3.4.

Linear Stability Analysis

23

The resultant expression (Eq.(37)) for K3 obtained from the CFD analysis has been inputted into the stability analysis to determine f1, (f1)D and (f1)N as a function of gas superficial velocity. The roots of these equations—which are meant to be the critical gas volume fraction(s), (G)crit, at which the system transitions from homogenous (stable) to heterogeneous (unstable) conditions—were then computed. Note that for a given constant superficial liquid velocity, UL = -0.134 m/s in this case, there will be a particular, (G)crit value corresponding to each superficial gas velocity, VL, value. This is because f1, (f1)D and (f1)N are all functions of VG; and also rB which change with VG as well. Hence, for a system with a fixed UL a curve of (G)crit versus VG will be generated. At low gas velocities this curve will normally be above the measured gas volume fraction, (G)meas, curve indicating that the system is operating under a stable (homogeneous) condition. At higher gas velocities the two curves will move closer together, and eventually at some superficial gas velocity value the two curves intersect and crossover and the experimental system will become unstable (heterogeneous). The gas volume fraction at which the intersection occurs will be the (G)crit condition for the given superficial liquid velocity, UL, that has remained constant. The (G)crit and (G)meas versus VG behaviour just described for a fixed UL is given in Fig. 16. However, before commenting on the final result the methodology used for determining (G)crit is firstly described with the assistance of Figs. 14 and 15, which have been produced using UL = -0.134 m/s, VG = -0.031 m/s, rB = 157 µm (from Fig. 5) and K3 = 3.3 (from Eq. (37)). Firstly, f1 (right side axis), and (f1)D and (f1)N (left side axis), have been plotted against G in Fig. 14. It can be clearly seen that all three functions have zero value (roots) at G between 0.3-0.4. Specifically, as shown in the expanded insert, the root for both f1 and (f1)N coincide with each other at G equal to 0.313; whereas the root for (f1)D occurs at G equal to 0.320. Whilst the difference is only about 2.2 percent, the difference of their influence on the physical system is significant. At G equal to 0.313 by virtue of (f1)N being equal to zero then f1 will have the same value; and the system will be at the point of becoming unstable. Possibly, the conditions within the system could change, i.e.

24

coalescence taking place, which could return to the system to a homogeneous (stable) condition operating at a new bubble size. Conversely, the root for (f1)D occurring at G equal to 0.320 represents a discontinuity in the f1 function, which presumably is no longer applicable beyond this value since the assumptions of 1D flow, constant bubble size, etc., are invalid. On this basis, a (G)crit value of 0.320 has been assumed for VG = -0.031 m/s; noting that the difference in taking the actual root of f1 (or(f1)N) is only around 2 percent. The close proximity of two roots for the stability parameter, f1, has been reported previously by Ghatage et al. (2014), with the comment that the first root (lowest G value) marked the transition from homogeneous to heterogeneous conditions. They also state that the second root, which should result in the transition back to the homogeneous condition, doesn’t actually occur in the real system since the phenomena associated with the heterogeneous conditions, i.e. coalescence, flow recirculation, etc., are irreversible. The results from this study are consistent with these observations, except that the second root is the consequence of a discontinuity of f1 due to a zero value in the denominator function, (f1)D. The hypothesis, therefore is that for the results presented in this study that the (G)crit value of 0.313, resulting from the root for f1—or(f1)N to be precise—marks the actual transition from homogeneous to heterogeneous conditions, but it is the root at (G)crit value of 0.320 that marks the limiting condition beyond which the basic assumptions that underpin the linear stability analysis no longer apply. As discussed in the modelling approach outlined in Section 2, there exist a number of other roots for the f1, (f1)N and (f1)D functions, and it’s worth considering their meaning in terms of the homogeneous-to-heterogeneous transition of the system. To assist in the discussion the (f1)N and (f1)D graphs given in Fig. 14 have been replotted in Fig. 15 for G in the range -0.2≤G≤1.2. The reason for doing this is to highlight the different types of roots obtained for both (f1)N and (f1)D, namely those that cross (open circles) and those that just touch (close) the horizontal axis. Given that the polynomial function for (f1)D is of order 20, there will almost certainly be additional roots outside the

25

range of G values shown in Fig. 15. However, they will have no relevance to real gas-liquid systems where the gas volume fraction is bounded between 0 and 1. For this reason the functional behaviour of (f1)N and (f1)D beyond these limits has been excluded from the following discussion. It can be seen that there are roots at G equal to 0 and 1, and are generated by the (G)p and (1-G)q terms given in Eqs. (14) and (19). These roots are at both limits of real values that can be possibly attained for the gas volume fraction, and it’s not surprising that the mathematical formulation of the stability problem exhibits zero values in both the numerator and denominator functions. The second touching root shown in Fig. 15 is at G equal to 0.188 which is equivalent to the no slip condition between the gas and liquid phases, and specifically is generated by the [VG(1G)–ULG]r terms given in Eqs. (14) and (19). The implication of this result is that the system will become unstable once the gas volume of the system reaches the no slip value of VG/(VG+UL). If so, then as shown in Fig. 3 the system should be unstable for all values of VG since the measured volume fraction curve is always above the corresponding no slip curve—which is what you would expect for a downflowing gas-liquid system where the buoyancy of the bubbles would always create a slip condition and result in a higher gas volume fraction value. The experimental observation, however, is that the system can operate in the bubbly (homogeneous) flow regime beyond the critical no slip volume fraction condition. Whilst there might actually be an instability that causes coalescences the system is able to re-establish a stable condition corresponding to a new equilibrium (larger) bubble size. Evidence of a gradual increase in bubble size, in response to overcoming the critical no slip volume fraction condition for each VG value, is shown in Fig. 3. Consequently, more gas can be inputted into the system but only to the point that any instability can be overcome a gradual increase in bubble size. In Fig. 15 it can be seen that the next root (instability) for (f1)N and (f1)D actually crosses the horizontal axis at (G)crit equal to 0.313 and 0.320, respectively. The actual meaning of these two values—especially being so close together—has been discussed previously, with the implied

26

assumption that the root for (f1)D being a good measure of the actual value (G)crit. The mathematical origin of this instability point is the polynomial expression within Eq. (19), which is generated by the X2-5 terms given by Eqs. (20)-(23), respectively. These terms comprise the variables K0-3, dB and CV0 that provide specific information about the operational conditions of the system, i.e. fluctuation velocity of the incoming liquid through K3. Therefore, it is thought that the numerical value of (G)crit generated by Eq. (19) from the X2-5 terms is more likely to be the actual critical gas volume fraction. It is important to note that the analysis is based on the assumption that the instability created by the no slip condition can be firstly overcome by balancing coalescence and breakup phenomena to reaching a new equilibrium bubble size. The stability function terms presented in Figs. 14 and 15 are for a given liquid superficial velocity and a single gas superficial velocity. The outcome from those figures is an (G)crit value that corresponds to that particular VG value. Hence, if (G)meas is less than (G)crit then the system will be stable at that gas velocity. The methodology was repeated for each of the experimental VG values, and the resultant (G)meas and (G)crit curves are presented in Fig. 16; noting that two (G)crit curves are given for (1) K3 = 0 (▪▪▪▪), and (2) K3 (▬) given by Eq. (34). The reason for plotting the critical gas volume fraction graphs for the two K3 values is the assess the influence if velocity fluctuations of the inlet liquid stream coming from the Mixing Zone and into the Two Phase Pipe Flow Zone where the stability measurements were made. It can be seen that for K3 = 0 (▪▪▪▪), which gives the condition that there is no inlet liquid velocity fluctuation, the (G)crit curve lies above that for (G)meas over the entire VG range. On this basis, the system will remain in the homogeneous condition without any transition to heterogeneous flow taking place. This prediction is not supported by visual observation, where churn-turbulent-like flow was observed for G above about 0.3. The second (G)crit curve (▬) is where K3 has been allowed to vary in accordance with Eq. (34), which is based on the CFD analysis. It can be seen that the (G)crit curve falls above the

27

(G)meas curve, indicating stable or homogeneous flow, until there is a cross-over that marks the onset of unstable or heterogeneous flow. The cross-over takes place at VG = 0.079 m/s, which corresponds approximately with the VG = 0.075 m/s shown in Fig. 5 where there is a change in slope of the computed bubble size. The conclusion, therefore, is that the stability analysis can be applied to the Pipe Flow Zone of a plunging liquid jet bubble column provided the inlet liquid velocity fluctuation is properly considered.

4.

Conclusions Linear stability analysis has been applied to the Two Phase Pipe Flow Zone in a plunging

liquid jet bubble column. In order to do so, drift-flux analysis has been applied to the experimental data to obtain an estimate of the bubble size as a function of gas and liquid superficial velocities; whilst CFD has been used to quantify the fluctuating velocity of the liquid exiting the Mixing Zone. The gas volume fraction versus gas superficial velocity curve showed no obvious sign of transition from homogeneous-to-heterogeneous flow even at gas volume fractions above 0.4. However, the complementary graph of bubble size exhibited a change in slope at VG  0.075 m/s, possibly due to a transition in flow condition. The linear stability analysis gave a corresponding VG value of a 0.079 m/s, based on the selection of a critical volume fraction corresponding to the discontinuity of the linear stability function, f1. In applying the analysis it was assumed that instability of the two phase flow did not take place at no slip flow condition, i.e. VG/(VG+UL). It was also found that the instability of the system was influenced by the velocity fluctuation of the inlet liquid stream, which had not been considered in previous analysis.

References

28

Anderson, T.B., Jackson, R., 1969. A fluid mechanical description of fluidized beds – comparison of theory and experiment. Industrial & Engineering Chemistry Fundamentals 8, 137-144. Batchelor, G. K., 1988. A new theory of the instability of a uniform fluidized bed. Journal of Fluid Mechanics 193, 75-110. Bhole, M.R., Joshi, J.B., 2005. Stability analysis of bubble columns: predictions for regime transition. Chemical Engineering Science 60, 4493-4507. Clift, R., Grace, J.R., Webber, M.E., 1978. Bubbles, droplets and particles, Academic Press. Ekambara, K., Dhotre, M.T., Joshi, J.B., 2005. CFD simulation of bubble column reactors:1D, 2D and 3D approach. Chemical Engineering Science 60, 6733-6746. Evans, G.M., 1990. A study of plunging jet column. PhD dissertation, University of Newcastle, Australia. Ghatage, S.V., Bhole, M.R., Padhiyar, N., Joshi, J.B., Evans, G.M., 2014. Prediction of regime transition in three-phase sparged reactors using linear stability analysis. Chemical Engineering Journal 235, 307-330. Gibilaro, L.G., Hossain, I., Foscolo, P.U., 1986. Aggregate behaviour of liquid fluidized beds. Canadian Journal of Chemical Engineering 64, 931-938. Jackson, R., 1963. The mechanics of fluidized beds: I. The stability of the state of uniform fluidization, Transactions of Institution of Chemical Engineering 41, 13-21. Joshi, J.B., Dinkar, M., Deshpande, N.S., Phanikumar, D.V., 2001. Hydrodynamic stability of multiphase reactors. Advances in Chemical Engineering 26, 1-130. Kynch, G.J., 1952. A theory of sedimentation. Transactions of Faraday Society 48, 166-176. Lahey, R.T.Jr., 1991. Void wave propagation phenomena in two-phase flow. AIChE Journal 37, 123-135. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Computational Methods in Applied Mechanics and Engineering 3, 269-289.

29

Monahan, S.M., Vitankar, V.S., Fox, R.O., 2005. CFD predictions for flow regime transitions in bubble columns. AIChE Journal 51, 1897–1923. Nedeltchev, S., Shaikh, A., 2013. A new method for identification of the main transition velocities in multiphase reactors based on information entropy theory. Chemical Engineering Science 100, 2-14. Upadhyay, R.K., Kaim, J., Roy, S., 2009. Investigation of downflow bubble columns: Experiments and modeling. Journal of Chemical Engineering of Japan 42, 156-161. Shaikh, A., Al-Dahhan, M.H., 2007. A review on flow regime transition in bubble columns, International Journal Chemical Reactor Engineering 5, 1-68. Sundaresan, S., 2003. Instabilities in fluidized beds. Annual Reviews of Fluid Mechanics 35, 63-88. Wallis, G.B., 1962. A simplified one-dimensional representation of two-component vertical flow and its application to batch sedimentation, Symposium on the Interaction between Fluids and Particles, London, 9-16.

Fig. 1

Typical stability plot

Fig. 2

Schematic of liquid velocity fluctuation around bubble showing K3

Fig. 3

(G)meas (○) and (G)no slip (----) versus VG from Evans (1990) [UL = -0.134 m/s, D = 0.044 m, DN = 0.0047 m, Vj = 11.5 m/s, air-water]

Fig. 4

Drift-flux plot showing experimental () V’G versus and G [Curve for UL = -0.134 m/s, VG = -0.031 m/s, rB = 157 µm]

Fig. 5

Experimental (○) G and calculated () rB versus VG [UL = -0.134 m/s]

Fig. 6

Computational geometry for plunging liquid jet column

Fig. 7

Computed G contours: VG = (A) -0.028, (B) -0.039, (B) -0.069 m/s [UL = -0.134 m/s, rB values taken from Fig. 5]

Fig. 8

Experimental () and computed average G versus z [UL = -0.134 m/s, experimental G at z = 0.950-1.000 m]

30

Fig. 9

Computed average k versus z [UL = -0.134 m/s]

Fig. 10

Computed average  versus z [UL = -0.134 m/s]

Fig. 11

Computed average K3 versus z [UL = -0.134 m/s]

Fig. 12

Computed average K3 at Mixing Zone exit and fitted curve (eq. 37) versus VG [UL = 0.134 m/s]

Fig. 13

Computed K3 at Mixing Zone exit (eq. 37) versus VG with cut-off set at 50 [Open circles (○) correspond to data points given in Fig. 3]

Fig. 14

f1, (f1)D and (f1)N versus G [UL = -0.134 m/s, VG = -0.031 m/s, rB = 157 µm, K3 = 3.3]

Fig. 15

(f1)D and (f1)N versus G [From Eq. 10, with UL = -0.134 m/s, VG = -0.031 m/s, rB = 157 µm, K3 = 3.3]

Fig. 16

Measured (○) G and computed (G)crit versus VG [UL = -0.134 m/s, (▬) K3 = Eq. (34), (▪▪▪▪) K3 = 0]

31

f1

1

2

2.0

stable

1.5

0.5 0.0 0.25 -0.5 -1.0

Fig. 1

Typical stability plot

0.30

ϵG

0.35

0.40

unstable

f1

1.0

32

v’

bubble

(u’)inlet=K3v’

virtual mass bulk

Fig. 2

Schematic of liquid velocity fluctuation around bubble showing K3

33

ϵG

no slip

0.6 0.5

G

0.4

0.3 0.2 0.1 0.0 0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.16

VG (m/s) Fig. 3

(G)meas (○) and (G)no slip (----) versus VG from Evans (1990) [UL = -0.134 m/s, D = 0.044 m, DN = 0.0047 m, Vj = 11.5 m/s, air-water]

34

Drift flux curve

0.018

measurement

intersection

0.016

V'G (m/s)

0.014 0.012 0.010 0.008 0.006 0.004 0.002 0.000 0.0

0.2

0.4

ϵG

0.6

0.8

1.0

Fig. 4 Fig. 4

Drift-flux plot showing experimental () V’G versus and G [Curve for UL = -0.134 m/s, VG = -0.031 m/s, rB = 157 µm]

35

ϵG

0.6

rB

0.0008 0.0007

0.5

0.0006 0.0005

0.3

0.0004 0.0003

0.2

0.0002 0.1

0.0001

0.0 0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.0000 0.16

VG (m/s) Fig. 5 Fig. 5

Experimental (○) G and calculated () rB versus VG [UL = -0.134 m/s]

rb (m)

ϵG

0.4

36

Gas IN

Liquid IN

Wall

Wall

A

A

OUTLET Fig. 6

Computational geometry for plunging liquid jet column

37

(A)

Fig. 7

(B)

(C)

Computed G contours: VG = (A) -0.028, (B) -0.039, (B) -0.069 m/s [UL = -0.134 m/s, rB values taken from Fig. 5]

38

0.028 m/s

0.049 m/s

0.069 m/s

(ϵG)exp

1.0

ϵG (-)

0.8

0.6

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

z (m)

Fig. 8

Experimental () and computed average G versus z [UL = -0.134 m/s, experimental G at z = 0.950-1.000 m]

39

0.028 m/s

0.049 m/s

0.069 m/s

2.0

k (m2/s2)

1.5

1.0

0.5

0.0 0.0

0.2

0.4

0.6

z (m)

Fig. 9

Computed average k versus z [UL = -0.134 m/s]

0.8

1.0

40

0.028 m/s

0.049 m/s

0.069 m/s

0.069

0.049

0.028

1000.00

ε (m2/s3)

100.00

10.00

1.00

0.10

0.01 0.0

0.2

0.4

0.6

z (m) Fig. 10

Computed average  versus z [UL = -0.134 m/s]

0.8

1.0

41

0.028 m/s

0.049 m/s

0.069 m/s

24 22 20 18

K3 (-)

16 14 12 10 8 6 4 2 0 0.0

0.2

0.4

0.6

z (m) Fig. 11

Fig. 11

Computed average K3 versus z [UL = -0.134 m/s]

0.8

1.0

42

(K3)MZexit

fitted

extrapolated

16 14

(K3)MZexit (-)

12 10 8 6 4 2 0 0.000

0.020

0.040

0.060

0.080

VG (m/s)

Fig. 12

Computed average K3 at Mixing Zone exit and fitted curve (eq. 37) versus VG [UL = 0.134 m/s]

43

60 cut-off value 50

K3

40 30

20 10 0 0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.16

VG(m/s) Fig. 13 Fig. 13

Computed K3 at Mixing Zone exit (Eq. 37) versus VG with cut-off set at 50 [Open circles (○) correspond to data points given in Fig. 3]

44

f1(D)

f1(N)

f1

f1 (discontinuity) 8

0.E+00 0.0

0.2

0.4

0.6

ϵG

-1.E+17

0.8

1.0

6 4

f1(N)

f1

f1 (discontinuity)

08

8.E+15

6

6.E+15

-3.E+17 -4.E+17

4.E+15

4 -2

2.E+15

2

0.E+00 0.25 -2.E+15 -4.E+15

-6.E+15

-5.E+17

Fig. 14

f1

1.E+16

f1 zero

ϵG

0.30

-4

0 0.35 -2

f1

-2.E+17

(f1)D, (f1)N

(f1)D, (f1)N

2 f1(D)

-4 -6

-8.E+15

-6

-1.E+16

-8 -8

f1, (f1)D and (f1)N versus G [UL = -0.134 m/s, VG = -0.031 m/s, rB = 157 µm, K3 = 3.3]

45

f1(D)

crossesf1(D) axis 1.E+16

5.E+10

no slip touches F1(D) root axis

f1(N)f1(N) F1(N) root

8.E+15 6.E+15 4.E+15

3.E+10

(f1)D, (f1)N

(f1)D

2.E+15

1.E+10

0.E+00 0.05 -2.E+15

0.10

0.15

0.4 ϵ G

0.6

0.20

0.25

0.30

0.35

-4.E+15

-6.E+15

-1.E+10

-0.2

0.0

-8.E+15

0.2 -1.E+16

0.8

1.0

1.2

-3.E+10

-5.E+10

Fig. 15

(f1)D and (f1)N versus G [From Eq. 10, with UL = -0.134 m/s, VG = -0.031 m/s, rB = 157 µm, K3 = 3.3]

46

ϵG

(ϵG)crit

instability

K3=0

0.8 0.7

G, (G)crit

0.6 0.5 0.4 0.3 0.2 0.1 0.0 0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.16

VG (m/s) Fig. 16

Measured (○) G and computed (G)crit curve versus VG [UL = -0.134 m/s, (▬) K3 = Eq. (34), (▪▪▪▪) K3 = 0]

47

f1(D)

f1(N)

f1

f1 (discontinuity)

f1 zero

1.E+16

8

8.E+15

6

6.E+15

4 4.E+15

(f1)D, (f1)N

0.E+00 0.25 -2.E+15

ϵG

0.30

0 0.35 -2

-4.E+15 -4

-6.E+15 -8.E+15

-6

-1.E+16

-8

Fig. 14

Inserted graph

Highlights 

Analyzed instabilities occurring in plunging jet bubble column.



CFD simulations predicted the phase fluctuations.



Applied stability criteria based on the linear stability analysis.



Significance of each of the roots that constitute the stability criteria.



Good agreement with published experimental results.

f1

2

2.E+15