Detailed Investigations of the Countercurrent Multiphase (Gas–Liquid

Apr 22, 2014 - Logistics Engineering College, Shanghai Maritime University, Shanghai 200135, China. ‡. Department of Process Dynamics and Operation,...
0 downloads 0 Views 5MB Size
Article pubs.acs.org/IECR

Detailed Investigations of the Countercurrent Multiphase (Gas− Liquid and Gas−Liquid−Liquid) Flow Behavior by Three-Dimensional Computational Fluid Dynamics Simulations YuanYuan Xu,*,† Ming Zhao,† Steve Paschke,‡ and Guenter Wozny‡ †

Logistics Engineering College, Shanghai Maritime University, Shanghai 200135, China Department of Process Dynamics and Operation, TU Berlin, 10623 Berlin, Germany



S Supporting Information *

ABSTRACT: This paper focuses on the detailed investigations of the (multi-) liquid flow behaviors on an inclined stainless steel plate in simulations and experiments. Simulations are carried out with FLUENT 6.3; a three-dimensional computational fluid dynamics model considering the gravity, surface tension, and local drag force is developed and presented. Experimental data is obtained with optical methods. One main aspect of the paper’s presentation is the validation of the developed model for the case of one-liquid phase with a countercurrent gas phase, where the parameters (e.g., the thickness and the velocity profile of the liquid phase) are compared with the simulation results and the experimental data. The other key part of the presentation is the numerical investigations with the developed model for the case of two immiscible liquid phases with a countercurrent gas phase. Results show that with the current gas flow rate, the complex flow behavior strongly depends on the (multi-) liquid flow rates and the feed sequence. Furthermore, there is a high stability between the two liquids. All of these data lead to a better understanding of the complex fluid dynamics in three-phase distillation processes. Because the change of the liquid flow behavior takes place only in a few inches, numerous investigations were first carried out on a single vertical plate with only one liquid phase in experiments.7−15 At this stage, some qualitative conclusions were drawn from the experimental observation and various correlations were analyzed with experimental data. To anticipate future work studying the liquid flow behavior inside packed columns, the experimental approach could not be followed because it is almost impossible to look into the packings. In recent years, with the development of computer technology, Computational Fluid Dynamics (CFD) studies on the local flow behaviors of liquid phases have become more and more popular.16−24 Nevertheless, in most of the literature, CFD simulations were carried out by simplifying the problem into a two-dimensional model18,19 and focusing on one-liquid phase.20−22 Furthermore, the interactions between the liquid and gas phases were neglected and therefore there were rare countercurrent simulations in the literature.23−25 However, in real industrial operations, the gas is in countercurrent flow to the liquid phases. The interaction between the gas and liquid phases is the key to successfully describing the flow behavior because it could lead to different flow patterns, such as reverse flow or entrainment. Therefore, the interaction between phases should be considered in the model. In addition, the local flow behavior such as film breakup with droplets and rivulets flow exhibits a highly three-dimensional characteristic; a twodimensional simulation could not capture the real information

1. INTRODUCTION During the last few decades, structured packed columns have been widely used in distillations or absorptions. This is due to their higher loading capacities and separation efficiencies as well as the lower pressure drop compared to traditional tray columns. However, these advantages are widely found in twophase distillations or absorptions. In the situation in which there is more than one liquid phase within distillation columns, the separation efficiencies become uncertain. This is mainly related to the fact that the liquid flow behavior changes dramatically with the appearance of a second liquid phase, which leads to the variations of the interfacial area for mass transfer, so that the separation efficiency changes greatly as well. This shows that the liquid flow behavior is an important factor affecting the separation efficiencies. In addition, it seems reasonable to consider the local flow behavior into the column models.1,2 To develop such a model, systematic investigations on twophase and three-phase packed distillations have been carried out from both the microscopic view (to study the liquid flow behavior on packings) and the macroscopic view (to make experiments on a laboratory-scale packed column) at the Institute of Process Engineering, TU Berlin. Chen et al.3,4 initially developed a three-phase mass-transfer model for packed distillations combined with the fluid dynamic parameters and obtained some valuable results, which is really an encouragement to further investigate the liquid flow behavior on the packing surface. As pointed out by Taylor,5 the work by Repke et al.1,2,6 is of particular interest because their results, such as film breakup and rivulet flow, have been used as fundamental assumptions in developing a mass-transfer model. © 2014 American Chemical Society

Received: Revised: Accepted: Published: 7797

January 4, 2014 March 19, 2014 April 4, 2014 April 22, 2014 dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

In each computational cell, the volume fractions of all phases sum to unity. Therefore, for the pth phase, its volume fraction γp has three situations: (a) γp = 0: the current cell is empty for the pth phase; (b) γp = 1: the current cell is full of the pth phase; (c) 0 < γp < 1: the current cell is the interface between the pth phase and one or more other phases. 2.1. CFD Equations. In general, the fluid with the lightest density is denoted as the primary phase in the VOF model; the other fluids are the secondary phases. Therefore, in our cases, gas is the primary phase and liquids are the secondary phases. The tracking of the interface between the phases is accomplished by the solution of a continuity equation for the volume fraction of the liquid phases, which is

(see Figure 1); the effective interfacial area between the gas and liquid phases estimated by a two-dimensional model is

∂γq

Figure 1. Possible flow behavior in packings: (a) two-liquid phase flow in 2D, (b) two-liquid phase flow in 3D (one flows as film and the other flows as rivulets and droplets).

∂t

+ ∇·(γqv ⃗) = 0

q = L1, L2

(1)

while the volume fraction of the gas phase will be computed based on the constraint

unreliable.3,19 Besides, the appearance of the second liquid sometimes is inevitable in the process industry, e.g., in bioethanol and biodiesel columns. In all, a three-dimensional countercurrent CFD model for multiliquid phases is necessary, which is the objective of this paper. This work focuses on the development and validation of models provided by CFD software to investigate the local flow behavior of the liquid phase and to examine the flow parameters such as film thickness, specific mass-transfer area, film velocity, and structure. First, a three-dimensional transient CFD model which is based on the volume of fluid (VOF) method is developed. This model considers the gravity, surface tension, and drag force between a liquid phase and a countercurrent gas phase. With the developed model, the flow characteristics of one liquid phase flowing down an inclined and flat plate in the case without and with a countercurrent gas flow are investigated. On one hand, the local and time-dependent hydrodynamics of the liquid phase, such as the film thickness, surface velocity, and velocity profile of the liquid phase are analyzed and compared with the experimental data to validate the developed model. On the other hand, the influence of the countercurrent gas flow on the liquid phase is investigated. Second, the flow behavior of two immiscible liquid phases flowing down the same plate in the case without and with a countercurrent gas flow is investigated in detail, where liquid flow behavior and flow characteristics as well as the influence of the countercurrent gas phase on liquid phases are discussed qualitatively.

γG + γL1 + γL2 = 1

(2)

The momentum equations along the x-, y-, and z-direction can be generally expressed as ∂ (ρ v )⃗ + ∇· (ρm vv⃗ ⃗ ) = −∇P + ∇· [μm (∇v ⃗ + ∇v ⃗ T)] + F ⃗ ∂t m (3)

where ρm and μm are the volume-averaged mixture density and viscosity, respectively, for all phases and computed in the manner below: ρm = γGρG + γL1ρL1 + γL2ρL2 (4) μm = γGμG + γL1μL1 + γL2μL2

(5)

F⃗ in eq 3 is the momentum source term, which denotes the influence of the external forces such as the gravity and the surface tension as well as the drag force between a liquid and a countercurrent gas phases. It is worth mentioning that the situation simulated here is a fully developed flow; the two liquids flow in the same direction and a velocity distribution is adopted at the boundary condition to avoid the velocity gradient; and consequently, the drag force between the two liquids is neglected in the current stage. The gravity source term (ρmg)⃗ can be included into the equation directly as the gravity is a volume force; whereas the surface tension source term (F⃗σ) is used by the continuum surface model (CSF),26 where an interfacial delta function is adopted to change the surface tension from a surface force to a volume force. F⃗σ in a computational cell has the following form:

2. CFD MODELING The purpose of CFD simulation is to investigate the local-flow behavior of the liquid phase under different liquid and gas loads. Because the gas and the two immiscible liquids cannot interpenetrate, the position of the interface and the interaction between the fluids are of interest. Therefore, the VOF model with the geometric reconstruction scheme is adopted to resolve the movement of the fluid. The VOF model is a surfacetracking technique applied to a fixed Eulerian mesh. In this model, a single set of momentum equations is shared by the fluids and the resulting velocity field is shared among the phases. For each phase in the model, a variable (the volume fraction of the phase in a computational cell, γ) is introduced, which is tracked throughout the whole computational domain.

Fσ⃗ pq =



σpq

pairs pq , p < q

γpρp κq∇γq + γqρq κp∇γp 1 (ρ 2 p

+ ρq )

p , q = G, L1, L2 (6)

For a film flow with a countercurrent gas phase, the driving force for the liquid phase is gravity, whereas that for the gas phase is the pressure drop. Therefore, a friction pressure drop model is considered here to describe the drag force source term F⃗f between a liquid and the countercurrent gas phases.27 F⃗f in a computational cell is Ff⃗ =

∑ p = G; q = L1,L2

7798

( −apqf pq ρp (up̅ − uq i)|up̅ − uq i|)

(7)

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

Figure 2. Geometry, mesh, and boundary conditions adopted into simulations for different cases.

Figure 3. Velocity profiles at the inlet boundaries for ReL=f(u, δ), with the solid line representing higher liquid loads and the dotted line representing lower liquid loads: (a) one-liquid flow, (b) two-liquid flow.

⎞1/3 ⎛ 3·(μ /ρ )2 q q ⎜ δq = ⎜ ·Req⎟⎟ g sin · α ⎠ ⎝

where apq is the effective interfacial area per unit volume between a gas and a liquid phase, which can be calculated by apq =

2|∇γp||∇γq| |∇γp| + |∇γq|

(8)

f pq = 0.079Rep−0.25(1 + 115δq*N )

δq* =

Bo =

N = 3.95/(1.8 + 3.0/Bo)

(10)

where the Reynolds number Re, film thickness δ, and Bond number Bo are calculated by ρp up̅ dh μp Vq̇ b·μq /ρq

=

ρq uq̅ δq μq

δq [σpq/(ρq − ρp )g ]1/2

q = L1 or L2 (13)

dh [σpq/(ρq − ρp )g ]1/2

(14)

2.2. Simulation Strategies. It has been known that the local flow behavior of liquid phases including the film flow and film breakup with rivulet and droplet flow takes place in the range of only a few inches on packings. Therefore, in the first step of investigating fundamental details, the structured packing is substituted and simplified by a small flat smooth plate. For the case of one-liquid flow, the geometry for simulations resembles the experimental setup with the flow field of 60 × 50 mm2 and a 10 mm height of the channel. For the case of twoliquid flow, because of the increase of the liquid thickness (two liquids), the height of the channel for simulations is increased to 15 mm with the aim of reducing the effect of the outflow field and the influence of the gas pressure on liquid phases. The flow field is reduced in half along the width of the plate, namely to 60 × 25 mm2 in order to save computational time. The

p = G, q = L1 or L2 (9)

Req =

(12)

p = G, q = L1 or L2

and f pq is the drag force factor between a gas and a liquid phase, which strongly relies on the geometry of the setup. Because in the current stage our experimental setup is similar to a rectangular channel, the drag force factor proposed by Stephan and Mayinger28 seems to be suitable in considering the distillation and absorption processes, which is described as

Rep =

q = L1 or L 2

p = G, q = L1 or L2 (11) 7799

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

Figure 4. Experimental setups for the case without the countercurrent gas flow,30,31 where (a) one liquid flows in the relatively microscopic experiment and (b) two liquids flow in the relatively macroscopic experiment.

residuals. When each monitored variable fluctuates around a constant, it can be regarded as a quasi-stable state. Then the gas phase is induced countercurrently to the liquid phase(s). The influence of the gas phase on the liquid phase(s) could be investigated by gradually increasing the gas flow rate.

commercial simulation package Fluent 6.3 (ANSYS Inc.) is used in the processes. Because the mesh and the boundary conditions are sensitive to the simulation results, considerable work has been done to determine them. The information about the mesh and boundary conditions for one-liquid flow could be found in the literature.24 Figure 2 shows the geometry, mesh and boundary conditions used in simulations for two-liquid flow. To consider the effect of the wall adhesion when a liquid is in contact with the wall boundaries and to capture the interfacial phenomena between phases, the meshes in the direction of the plate and on both sides as well as around the interface are very fine. In addition, considering the different situations for the case without or with the countercurrent gas flow, the boundary conditions vary from case to case (see the table in Figure 2). To avoid numerical divergence caused by the velocity difference between two phases, a velocity profile is used at the inlet boundary. The velocity-inlet boundaries include the velocity profiles and the thicknesses of liquid phases, which are obtained by the Nusselt film theory at a given Reynolds number of the liquid phase.29 Figure 3 displays the velocity profiles at the inlet boundary in the case without the countercurrent gas phase in this work. It can be seen that velocity profiles for one-liquid flow are based on Nusselt theory and that for two-liquid flow are formed on the basis of both Nusselt theory and the plug flow. In addition, both the maximum velocity and film thickness decrease with the decrease of the Re-number. In the countercurrent cases, the pressure-inlet boundaries are set according to the velocities of the phases. For example, a pressure profile based on the velocity profile of the liquid phase is used at the inlet boundary; the pressure at the outlet boundary is determined by the desired velocity of the gas phase ugd. It is worth noting that the film velocity profiles and the pressure-inlet boundaries together with the drag force source term are integrated into FLUENT by user defined functions (UDF). Because the Re-numbers of the liquid phase in the considered cases are lower than 250 and the velocities of the countercurrent gas phase are not very high (lower than 6 m/s), turbulence is not obvious; therefore, it could be treated as laminar flow. Simulations start with the case(s) without the countercurrent gas flow, where the mass flow rate for the whole domain and the force on the plate are monitored beside

3. EXPERIMENTAL SETUP To validate the developed CFD model, experiments are carried out in two ways according to the investigation area on the plate: the relatively macroscopic experiment and the relatively microscopic experiment. The relatively macroscopic experiment concentrates the research on the whole plate, where the global parameters such as the wetted area of the plate, liquid thickness, and liquid surface velocity along the flow direction are measured and analyzed;30 the relatively microscopic experiment focuses the research on a small region on the plate, where the local parameters such as the velocity profile along the liquid thickness are scanned and investigated layer-by-layer.31 Figure 4 shows the experimental setups for the case without the countercurrent gas flow. Figure 4a displays the microscopic experimental setup and Figure 4b shows the macroscopic experimental setup. In experiments, the plate is a stainless steel plate with a flow field of 80 × 50 mm2 and the inclination angle of 60°, which is a representative angle of structured packings (X Mellapak). The overflow weir or the feed tube delivers the liquid(s) onto the plate. In macroscopic experiments, a laserinduced fluorescence (LIF) method is adopted to measure the liquid thickness; the surface velocity along the flow direction is determined by the particle tracking velocity (PTV) method. In microscopic experiments, the particle image velocimetry (PIV) method is used from the top or the bottom of the plate to measure the velocity profile along the film thickness. Detailed operations for the experimental process can be found in the literature30,31. It is worth mentioning that experimental measurements mainly focus on the film flow because the optical measuring system is very sensitive to the curvature of the rivulet flow. 4. RESULTS Simulations with the developed model were carried out for different systems and different flow patterns. The physical properties of the material used in simulations are presented in Table 1, where the contact angle θ among three phases is very 7800

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

Table 1. Physical Properties of the Materials in Simulations (T = 25°C, P = 1 bar) ρ (kg m−3) air toluene water

ν (m2 s−1)

σL−G (N m−1)

σL−L (N m−1)

θL−G (deg)

θL−L (deg)

− 0.029 0.073

− 0.0361 −

− 82 702

− 1242 −

−5

1.545 × 10 6.82 × 10−7 8.9 × 10−7

1.185 867 997

sensitive to the simulation results32 and is defined as the angle between the wall and the tangent to the interface at the wall. The schematic diagram for measuring the contact angle is shown in Figure 5. For each pair of phases, the phase with lighter density is specified as the primary phase.33 The materials and the parameters for the simulations are listed in Table 2.

4.1.1. Specific Wetted Area. The investigations of the specific wetted area are carried out for one-liquid flow without the countercurrent gas phase. Figure 6 shows the comparison of

Figure 5. Schematic diagram for measuring the contact angle between fluid−fluid phases.

Figure 6. Comparison of the specific wetted areas between experimental data in ref 1 and simulation results for case 1 and case 2.

Table 2. Cases for Simulations case

system

1

water, air

2

toluene, air water, air

3 4 5 6 7 8 9

ReL

ugd (m/s)

224, 168, 140, − 112, 84, 56 220, 147, 74, − 37 224 3, 5, 6

toluene, air water, air

147

3, 5

112

3, 5, 6

water, toluene, air water, toluene, air water, toluene, air water, toluene, air

ReW = 224, ReT =147



one liquid flowing on the plate one liquid flowing on the plate one liquid with countercurrent gas one liquid with countercurrent gas one liquid with countercurrent gas two liquids flowing on the plate (water is fed onto toluene)

ReW = 84, ReT − =165

two liquids flowing on the plate (toluene is fed onto water)

ReW = 224, ReT =147

two liquids with countercurrent gas (water is fed onto toluene)

3, 5

ReW = 84, ReT 3,5 =165

the specific wetted areas between experimental data and simulation results for case 1 and 2. For case 1, the specific wetted area is 1 at a higher ReL, which means water flows on the plate as a closed film; it falls below 1 when the ReL is lower than a certain value, which indicates that the water flow rate on the plate cannot support a closed film and the film breaks up. Furthermore, the specific wetted area decreases with the decrease of the ReL. These results show good agreement with the experiments as well. However, toluene flows always as a closed film under the current investigation conditions for case 2. Therefore, the flow behavior of the liquid phase strongly depends on the liquid flow rate and the properties of the liquid phase. Additionally, to show further the flow behavior of the liquid phase such as film flow or rivulet flow, a qualitative comparison between experimental data and simulation results is presented in Figure 7. These pictures display the same flow behavior between simulation and experiment, which also illustrates that the CFD simulation is able to predict the intricate phenomena on packings. 4.1.2. Liquid Thickness. The flow behavior of the liquid phase including film flow, rivulet flow, and droplet flow can be characterized by the liquid thickness and the velocity. In this work, studies are emphasized on the film flow and rivulet flow in the case without and with the countercurrent gas flow. Figure 8 shows the transient liquid film thickness over the plate width for case 1 with ReL = 224 (ZY plane, X = 50 mm). A good agreement between the measured profile by PIV and CFD simulations is shown in Figure 8. Furthermore, both simulation and experimental profiles display the formation of meniscuses on both sides. This is due to the support on both sides in experiments as well as the consideration of the contact angle at wall boundaries in simulations. Figure 9 displays the transient film thickness along the flow direction (XY plane, Z = 25 mm) against gas velocities for case 3 and case 4. It is worth noting that because of the limiting of boundary conditions, the pressure at the outlet boundary must

flow pattern

two liquids with countercurrent gas (toluene is fed onto water)

4.1. Two-Phase Flow (One-Liquid Phase Flow). First simulations were carried out for one liquid phase flowing on an inclined plate in the case without the countercurrent gas flow. Here film flow and film breakup with rivulets and droplets were investigated for case 1 and case 2. A parameter, named specific wetted area, which is defined as the ratio of the wetted plate area to the geometrical plate area, is introduced to study the flow behaviors. Validations between the simulation results and the experimental data were conducted as well. Afterward, simulations for film flow and rivulet flow in the case with the countercurrent gas flow are carried out separately. The influences of the countercurrent gas phase on the film flow and rivulet flow were investigated, and the liquid thickness and velocity are discussed in detail. 7801

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

Figure 7. Morphological descriptions of the flow behaviors for case 1: simulations (top) and experiments (bottom).

countercurrent gas flow, the film thickness is raised with the introduction of the countercurrent gas phase and increases with the growth of the gas velocity. This can be attributed to the drag force caused by the different velocities between the liquid and gas phase. The drag force tends to reduce the water flowing downward, which increases the liquid holdup on the plate so that the film thickness is raised. Figure 10 presents the formation process of the rivulet flow as well as the influence of the countercurrent gas phase on the

Figure 8. Comparison of the transient liquid thickness over the plate width between our simulations and the experiments in ref 1 for case 1 with ReL = 224 (ZY plane; X = 50 mm).

Figure 10. Instant rivulet thicknesses along flow direction (Z = 25 mm) against gas velocities for case 5.

rivulet flow. On one hand, when the film inlet condition such as the liquid loading is not able to maintain the film flow on the whole plate, the rivulet flow is formed. Meanwhile, a maximum hump appears on the plate because of the film breakup and the merging of the rims. On the other hand, similar to the film flow with the countercurrent gas phase, a part of the liquid on the maximum hump is prevented from flowing downward by the countercurrent drag force so that the liquid holdup on the plate increases; therefore, the height of maximum hump is raised. Simultaneously, the drag force directly carries the maximum hump upward, overcoming the gravity, until the forces among the gravity, surface tension, and drag force come to a new balance. Consequently, with the increase of air velocities, the maximum hump is both raised and drawn back. In order to further illustrate the variation of the liquid thickness and to validate this developed model with the experimental data, additional simulations corresponding to the experiment in ref 15 were carried out. In these cases, the liquid

Figure 9. Instant film thicknesses along flow direction (Z = 25 mm) against gas velocities. Lines with solid symbols represent case 3, and lines with hollow symbols represent case 4.

overcome the pressure generated by the liquid velocity and the liquid thickness. Therefore, the effective gas velocity is smaller than the desired gas velocity. Besides, with the same desired gas velocity ugd, slight differences of the effective gas velocity uge for different cases could be found as well. First, it can be seen in Figure 9 that in the case without the countercurrent gas flow (ugd = uge = 0 m/s), the waves on the film surface due to the gravity and the surface tension can be found in both cases. Furthermore, the fluctuation of the film flow for case 3 is higher than that for case 4, which indicates that the fluctuation of the film flow has a relationship with the material and the liquid flow rate. Second, it also can be seen that in the case with the 7802

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

value. This could be explained by the lower edge of the plate in experiments, which dams up the liquid there. It also can be seen in Figure 12 that the surface velocity in experiments obtained by the feed tube is higher than that obtained by the overflow weir. This is due to a smaller inlet region of the feed tube with the same flow rate. However, both the experiment and simulation results express the same conclusion: with the same flow rate, the rivulet flow runs faster than the film flow. The higher velocity means lower residence time, which is unfavorable to mass and heat transfer in distillations. To validate the model and to show the creditability of the simulation results, another simulation with a water−glycerol mixture and air on an inclined glass plate is carried out, which corresponds to the experiment in the ref 34. Figure 13 shows a

phase is water with surfactant and the gas phase is air. A parameter, named the mean film thickness, is calculated by z2

δ ̅ = (∑

x2

∑ δij)/((x2 − x1)(z 2 − z1))

i = z1 j = x1

(15)

A comparison of the mean film thickness against gas velocities between simulation and the experimental data is shown in Figure 11. Good agreement between them could be found.

Figure 11. Comparison of the mean film thickness against gas velocities between simulation and experimental data in ref 15.

Furthermore, both the experiment and simulation exhibit a similar tendency: the mean film thickness slightly increases with the growth of the gas velocities. It can be found that in simulations the effective gas velocity that triggers flooding is about 3.5 m/s, whereas in experiments it is about 3.31 m/s. 4.1.3. Velocity. Velocity is an important parameter for mass and heat transfer because it has a strong relationship with the residence time of the liquid phase in distillations. In this work, the investigations of the velocity are conducted in two ways: the macroscopic method and the microscopic method. Therefore, the surface velocity along the liquid flow direction and the velocity profiles along liquid thickness are investigated, respectively. Figure 12 shows the comparison of the surface velocity along liquid flow direction between experimental data (ReW ≈ 110)

Figure 13. Comparison of the velocity profiles along film thickness between simulations and the experimental data in ref 34 in the case without the countercurrent gas flow.

comparison of the velocity profile along liquid film thickness between simulation results and the experimental data in the case without countercurrent gas flow. It can be seen that simulation results and the experimental data are in good agreement. Because of the liquid viscous friction, the velocity of the liquid phase is very low near the plate and increases steadily with the distance from the plate. The profile meets a semiparabolic profile which has been theoretically analyzed by Nusselt.29 The influence of the countercurrent gas phase on the liquid phase is investigated in terms of the film flow and the rivulet flow. Figure 14 displays the local velocity profiles along film thickness against gas velocities for case 3 and case 4. In the case of uge = 0 m/s, both the velocity profiles in water film and in toluene film express semiparabolic distributions and the maximum velocities appear at the interface. When the air flows countercurrently to the liquid phase with a lower velocity (uge ≈ 2 m/s), drag forces caused by the different velocities between the gas and liquid phases prevent a part of the interfacial liquid from flowing downward, which slightly elevates the film thicknesses and reduces the liquid velocities. With the increase of air velocities, drag forces intensify as well. At a higher air velocity (uge > 4 m/s), higher drag forces further decrease the velocities at the liquid interface so that there are changes of the velocity profiles in the water and toluene films; the maximum velocities appear not at the interfaces but at a position under the liquid surface. Figure 15 presents the local velocity profiles of water rivulets flow against air velocities for case 5; Figure 15a shows the velocity in the x-direction, and Figure 15b shows the velocity in z-direction. As can be seen in Figure 15a, in the case of uge = 0 m/s, the velocity profile along the rivulet thickness expresses

Figure 12. Comparison of the surface velocity along the liquid flow direction between the experimental data (ReW ≈ 110) and CFD simulations (ReW = 112) for case 1.

and CFD simulations (ReW = 112) for case 1. It can be seen that the simulation results show good agreement with the experimental data obtained by the overflow weir, especially in the film flow region. The difference occurs at the end of the plate where the rivulet velocities in simulations still accelerate whereas those in the experiments seem to reach a constant 7803

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

Figure 14. Local velocity profiles along film thicknesses (X = 43 mm, Z = 25 mm) against gas velocities for case 3 (a) and case 4 (b).

Figure 15. Local velocity profiles along rivulet thickness (X = 43 mm, Z = 25 mm) against gas velocities for case 5. Panel a is the velocity in the xdirection, and panel b is the velocity in the z-direction.

Figure 16. Morphological descriptions of two-liquid flow with different flow rate and feed sequence.Panel a is for case 6 and panel b is for case 7.

7804

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

Figure 17. Liquid thicknesses over plate width at different x-positions against gas velocities. Panels a and b are for case 6 and case 8, and panels c and d are for case 7 and case 9.

4.2.1. Morphologic Description. The liquid flow rate and the feed sequence have impacts on the liquid flow behavior. Here the feed sequence is the order of the liquids that are delivered on the plate. Figure 16 displays the morphological descriptions of two-liquid flow with different flow rates and feed sequences for case 6 and case 7. It has been noted in the case of one-liquid flow that water or toluene flows as a closed film on the plate with a high flow rate of ReW = 224 or ReT = 147. However, the closed film flow behavior is disrupted when water is fed onto the toluene with the flow rates of ReW = 224 and ReT = 147: water flows as a nearly closed film, whereas toluene film breaks up with rivulet and flows below the water film (see Figure 16a). Figure 16b expresses the liquid flow behavior when one liquid flows with a high flow rate and the other one flows with a low flow rate. It has been shown in the case of one-liquid flow that water with the flow rate of ReW = 84 forms a rivulet flow. Here in the case of two-liquid flow, when toluene is fed onto water, water forms rivulets and droplets which are covered by or running on the toluene film. To sum up, Figure 17 indicates that the flow behavior of a two-liquid flow is not a superposition of two single-liquid flows. There is a strong interaction between the two liquids. 4.2.2. Liquid Thickness. Because the liquid flow behavior changes greatly on the plate in three-phase flow, the analysis of the liquid thickness is performed from the standpoint of different positions. Figure 17 describes the transient liquid thickness over the plate width at different x-positions for cases 6−9 (ZY plane). A higher gas velocity causes the breakup of the toluene film earlier (see Figure 17a). When water is fed onto the toluene phase with high flow rates, at the position of X = 32 mm, both liquid phases cover the whole plate width in the case without the countercurrent gas phase. However, at the same position, the toluene film breaks up with rivulet and flows

the structure: below a certain value the velocity profile is semiparabolic, whereas above the value the velocity is nearly a constant. The influence of the countercurrent air on the water rivulet flow is similar to that on the water film flow. At a lower velocity (uge = 2.1 m/s), the drag force prevents a part of water at the interface from flowing downward, which gives rise to an increase of the rivulet thickness accompanied by a decrease of the velocity in the rivulet. The drag force strengthens with the growth of air velocity. At a higher air velocity (uge = 4.5 m/s), the velocity profile becomes more complicated because of the drag force, gravity, and viscous force in the liquid phase. Under this circumstance, the velocity profile needs to be investigated in detail in the next step. Figure 15b indicates the z-velocity profile along the rivulet thickness, which exhibits the swing degree of the rivulet flow over the plate width. For example, in the case of uge = 3.7 m/s, the z-velocity along the rivulet thickness is obvious and unidirectional, which leads to the change of the axis and the swing of the rivulet flow. The Supporting Information contains a movie representing the morphological descriptions of the meandering rivulet flow against gas velocities. This swing of the rivulet flow illustrates again that a three-dimensional model is necessary to analyze the complicated flow behavior on packings. 4.2. Three-Phase Flow (Two-Liquid Phase Flow). In the process industries, the appearance of the second liquid sometimes is inevitable. The appearance of the second liquid makes the multiliquid flow behavior much more complex. On the basis of the studies on the flow behavior of one-liquid phase above, the investigations on multiliquid flow behavior of two immiscible liquids (water and toluene) over an inclined stainless steel plate are carried out. Similar to the two-phase (gas−liquid) flow, three-phase (gas−liquid−liquid) flow can also be characterized by the liquid thickness and the velocity. 7805

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

below the water film at uge = 4.0 m/s. Figure 17b displays that compared with the liquid thicknesses in the case without countercurrent gas flow, the rivulet thickness of the toluene phase is raised by the countercurrent gas phase at the position of X = 50 mm. This leads to an elevation of the water film flowing on the toluene phase. Panels c and d of Figure 17 illustrate that when toluene is fed onto the water phase, the water first forms a rivulet and flows below the toluene film on the upper part of the plate, and then the rivulet is entrained onto the toluene film on the lower part. Furthermore, the countercurrent gas phase increases the thickness of the water rivulet and the toluene film as well. 4.2.3. Velocity. As indicated above, the two liquids interact strongly with each other, which could be reflected by the surface velocity. Figure 18 shows the comparison between

Figure 19 expresses the average surface velocities along the flow direction against gas velocities for one-liquid and twoliquid flows, where the result for two-liquid flow comes from case 8 and that for one-liquid flow comes from ReW =224 or ReT =147. When water is fed onto the toluene phase, the average surface velocities for each liquid fall between the average surface velocities of each liquid with the same flow rate. Namely, the toluene flows faster in the two-liquid case than in the one-liquid case, whereas the water flows slower in the twoliquid case than that in the one-liquid case. Figure 19a,b also indicates that in the case of countercurrent gas flow, the change of the average surface velocities for toluene or water phases in the two-liquid case is smaller than that in the one-liquid case. In other words, the influence of the countercurrent gas flow on the liquids in the three-phase flow is not as obvious as that in the two-phase flow. These further illustrate that there is a strong interaction between the two liquids, which makes the twoliquid flow more stable. Similar to the studies on the liquid thickness, the analysis of the velocity profiles is conducted with regard to different positions. Figure 20 displays the local velocity profiles along liquid thicknesses against gas velocities at different positions for case 8 and case 9. When water is fed onto the toluene phase, the velocity profiles along liquid thicknesses resemble the film velocity profiles in the one-liquid case, which also express semiparabolic distributions. In the middle of the plate, the velocity profiles along toluene film thickness are almost unaffected by the countercurrent gas flow. The countercurrent gas just decreases the velocity profiles of the water phase (see Figure 20a). At the position near to the end of the plate, the velocity profiles for both liquids are reduced by the countercurrent drag force (see Figure 20b). Especially in the case of uge = 4.0 m/s, the thickness of the toluene rivulet is increased, which leads to a notable elevation of the water film thickness accompanied by the decrease of the velocity profiles along liquid thicknesses. When toluene is fed onto the water phase, the velocity profiles along liquid thicknesses become more complicated because of entrainment. In the middle of the plate, when the water rivulet flows below the toluene film, the velocity profiles along liquid thicknesses resemble the rivulet velocity profiles in the one-liquid case. The countercurrent gas phase decreases the velocities of the toluene film and has no obvious effect on the water phase (see Figure 20c). At the position near to the end of the plate, the water is entrained by the toluene phase and the two liquids blend into each other (see Figure 20d). In this case, the velocity profiles along the liquid thickness is a combination of the velocity profiles for film flow, rivulet

Figure 18. Comparison between experimental data and simulation results of the surface velocities for films flowing on an inclined plate without the countercurrent gas flow. Hollow symbols represent oneliquid flow, and solid symbols represent two-liquid flow.

experimental data and simulation results of the surface velocities for films flowing on an inclined plate in the case without the countercurrent gas flow. As can be seen in Figure 18, both the experimental data and simulation results express the same trend: surface velocities of the film flow in the twoliquid case are smaller than those in the one-liquid case. In simulations, it is also found that in the case of one-liquid phase flow, the surface velocities of the water film or toluene film change only in a very small range; whereas in the case of twoliquid phase flow, they have a wide range of variation. Especially in the case of two-liquid flow, the simulation results and the experimental data at ReT ≈ 160 show a good agreement.

Figure 19. Average surface velocities along flow direction against gas velocities. Hollow symbols represent one-liquid flow with ReW = 224 or ReT = 147; solid symbols represent two-liquid flow for case 8. 7806

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

Figure 20. Local velocity profiles along liquid thicknesses against gas velocities at different positions for different cases: (a) X = 30 mm, Z = 12.5 mm, case 8; (b) X = 55 mm, Z = 12.5 mm, case 8; (c) X = 30 mm, Z = 12.5 mm, case 9; (d) X = 55 mm, Z = 12.5 mm, case 9.

flow, and droplet flow, which need to be investigated and validated in future work.

On the basis of the studies of the flow behavior of one-liquid phase, investigations with the developed model on the multiliquid flow behavior for three-phase (gas−liquid−liquid) flow were performed. The flow behavior and characteristics of two immiscible liquids (water and toluene) as well as the influence of the countercurrent gas phase on liquid phases were analyzed and discussed. The results show that the flow behavior for the two immiscible liquids is very complicated: one forms film flow, whereas the other forms rivulets or droplets and flows below or above the film surface. Sometimes entrainment could be observed on the plate. The complex flow behavior strongly depends on the liquid flow rates and the feed sequence. Moreover, the two liquids strongly interact with each other, which leads to a high stability in the case of countercurrent gas flow. Therefore, the influence of the countercurrent gas flow on the liquids in three-phase flow is not as obvious as that in twophase flow. All the investigations in this work provide a fundamental knowledge of the flow performance on packings, which is essential for developing a predictive model in threephase distillations in the future.

5. CONCLUSIONS A three-dimensional CFD model has been developed for the investigations of the liquid flow behavior of one-liquid and twoliquid phases flow on inclined and flat plates in this work. This model, based on the VOF method, is developed by UDF to consider the gravity and the surface tension as well as the local drag force between a liquid and a countercurrent gas phases. With the developed model, detailed investigations on the liquid flow behavior for two-phase (gas−liquid) flow over an inclined stainless steel plate have been carried out for the first time. The possible flow behavior and characteristics for oneliquid phase such as film flow, film breakup with rivulet, and droplet flow as well as the influence of the countercurrent gas phase on the liquid phase are analyzed and discussed. The experimental data was obtained (1) from our own experiments in the case without the countercurrent gas flow and (2) from the literature in the case with the countercurrent gas flow. These comparisons were used to demonstrate the feasibility of the CFD simulations and to validate the developed model. Both simulation results and the experimental data indicate that the investigations on the liquid flow behavior have to be conducted in three dimensions and that the flow rate and the property of the liquid phases have important influences on the liquid flow behavior. Results also indicate that the countercurrent gas phase not only increases the liquid thickness and decreases the liquid velocity, including surface velocity and the velocity profiles, but also affects the flow manner of the rivulet flow so that a swing of the meandering rivulet could be observed on the plate.



ASSOCIATED CONTENT

S Supporting Information *

Movie representing the morphological descriptions of the meandering rivulet flow against gas velocities. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*: [email protected]. Notes

The authors declare no competing financial interest. 7807

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research



Article

(2) Hoffmann, A.; Ausner, I.; Repke, J.-U.; Wozny, G. Detailed investigation of multiphase (gas-liquid and gas-liquid-liquid) flow behavior on inclined plates. Chem. Eng. Res. Des. 2006, 84, 147. (3) Chen, L.; Repke, J.-U.; Wozny, G.; Wang, S. Q. Exploring the essence of three-phase packed distillation-substantial mass transfer computation. Ind. Eng. Chem. Res. 2010, 49, 822. (4) Chen, L.; Repke, J.-U.; Wozny, G.; Wang, S. Q. Extension of the mass transfer calculation for three-phase distillation in packed columnnonequilibrium model based parameter estimation. Ind. Eng. Chem. Res. 2009, 48, 7289. (5) Taylor, R. Distillation modeling after all these years: A view of the state of art-catching up with the past. Ind. Eng. Chem. Res. 2007, 46, 4349. (6) Repke, J.-U.; Ausner, I.; Paschke, S.; Hoffmann, A.; Wozny, G. On the Track to Understanding Three Phases in One Tower. Chem. Eng. Res. Des. 2007, 85, 50. (7) Karimi, G.; Kawaji, M. An Experimental study of freely falling films in a vertical tube. Chem. Eng. Sci. 1997, 53 (20), 3501−3512. (8) Moran, K.; Inumaru, J.; Kawaji, M. Instantaneous hydrodynamics of a laminar wave liquid film. Int. J. Multiphase Flow 1999, 28, 731. (9) Vlachogiannis, M.; Bontozoglou, V. Experiments on laminar film flow along a periodic wall. J. Fluid Mech. 2002, 457, 133. (10) Akio, M. Numerical simulation of wavy liquid film flowing down on a vertical wall and an inclined wall. Int. J. Therm Sci. 2000, 39, 1015. (11) Liu, J.; Schneider, J. B.; Gollub, J. P. Three-dimensional instabilities of film flows. Phys. Fluids 1995, 7, 55. (12) Zapke, A.; Kröger, D. G. Countercurrent gas−liquid flow in inclined and vertical ducts  I: Flow patterns, pressure drop characteristics and flooding. Int. J. Multiphase Flow 2000, 26, 1439. (13) Drosos, E. I. P.; Paras, S. V.; Karabelas, A. J. Counter-current gas-liquid flow in a vertical narrow channel-liquid film characteristics and flooding phenomena. Int. J. Multiphase Flow 2006, 32, 51. (14) Zhu, M.; Liu, C. J.; Zhang, W. W.; Yuan, X. G. Transport phenomena of falling liquid film flow on a plate with rectangular holes. Ind. Eng. Chem. Res. 2010, 49, 11724. (15) Roy, R. P.; Jain, S. A study of thin water film flow down an inclined plate without and with countercurrent air flow. Exp. Fluids 1989, 7, 318. (16) Tong, Z. Y.; Marek, A.; Hong, W. R.; Repke, J.-U. Experimental and numerical investigation on gravity-driven film flow over triangular corrugations. Ind. Eng. Chem. Res. 2013, 52, 15946. (17) Owens, S. A.; Perkins, M. R.; Eldridge, R. B. Computational fluid dynamics simulation of structured packing. Ind. Eng. Chem. Res. 2013, 52, 2032. (18) Gao, D.; Morley, N. B.; Dhir, V. Numerical simulation of wavy falling film flow using VOF method. J. Comput. Phys. 2003, 192, 624. (19) Gu, F.; Liu, C. J.; Yuan, X. G. CFD simulations of liquid film flow on inclined plates. Chem. Eng. Technol. 2004, 27, 1099. (20) Xu, Z. F.; Khoo, B. C.; Wijeysundera, N. E. Mass transfer across the falling film: Simulations and experiments. Chem. Eng. Sci. 2008, 63, 2559. (21) Ataki, A. Wetting of structured packing elements - CFD and experiment. Ph.D. Dissertation, TU Kaiserslautern, Kaiserslautern, Germany, 2006. (22) Szulczewska, B.; Zbicinski, I.; Gorak, A. Liquid flow on strucured packing: CFD simulations and experimental study. Chem. Eng. Technol. 2003, 26, 584. (23) Lakehal, D.; Fulgosi, M.; Yadigaroglu, G. Direct numerical simulation of turbulent heat transfer across a mobile, sheared gasliquid interface. J. Heat Trans. 2003, 125, 1129. (24) Xu, Y. Y.; Paschke, S.; Repke, J.-U.; Yuan, J.; Wozny, G. Portraying the Countercurrent Flow on Packings by Three-Dimensional Computational Fluid Dynamics Simulations. Chem. Eng. Technol. 2008, 31, 1445. (25) Sun, B.; He, L.; Liu, B. T.; Gu, F.; Liu, C. J. A new multi-scale model based on CFD and macroscopic calculation for corrugated structured packing column. AIChE J. 2013, 59, 3119. (26) Brackbill, J. U.; Kothe, D. B.; Zemach, C. A continuum method for modelling surface tension. J. Comput. Phys. 1992, 100, 335.

ACKNOWLEDGMENTS The authors acknowledge the Shanghai Natural Science Fund (Grant 13ZR1418900), Innovation Program of Shanghai Municipal Education Commission (Grant 14YZ105), National Natural Science Foundation of China (Grant 61302132), Specialized Research Fund for the Doctoral Program of Higher Education (Grant 2013312120001), and Science & Technology Program of Shanghai Maritime University (Grant 20130435) for financial support.



NOMENCLATURE a (m2/m3) = effective interfacial area per unit volume b (m) = width of the plate Bo (−) = bond number dh (m) = hydraulic diameter of the channel f (−) = drag force factor F⃗ (kg/m3 s2) = momentum source term g ⃗ (m/s2) = acceleration of gravity P (Pa), (bar) = pressure Re (−) = Reynolds number t (s) = time T (°C) = temperature u (m/s) = velocity in x-direction u̅ (m/s) = average velocity uge (m/s) = effective gas velocity ugd (m/s) = desired gas velocity v (m/s) = velocity in y-direction V̇ (m3/s) = flow rate v ⃗ (m/s) = velocity vector, v ⃗ = (u,v,w) x, X (mm) = coordinate in the direction of the liquid flow y, Y (mm) = coordinate in the direction of the liquid thickness z, Z (mm) = coordinate in the direction of the plate width α (deg) = inclination angle γ (−) = volume fraction δ (m) = film thickness δ̅ (m) = average film thickness δ* (-) = dimensionless ratio of film thickness to the Laplace length μ (kg m−1 s−1) = dynamic viscosity ρ (kg m−3) = density σ (N m−1) = surface tension θ (deg) = contact angle κ (−) = curvature of the interface between two phases ∇ = del operator, ∇ = ∂/∂x + ∂/∂y + ∂/∂z

Subscript and Superscript

G = gas phase L1 = first liquid phase L2 = second liquid phase m = mixture i = interfacial p, q = phase (gas, liquid 1, and liquid 2) f = drag force σ = surface tension W = water T = toluene



REFERENCES

(1) Hoffmann, A.; Ausner, I.; Repke, J.-U.; Wozny, G. Fluid dynamics in multiphase distillation processes in packed towers. Comput. Chem. Eng. 2005, 29, 1433. 7808

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809

Industrial & Engineering Chemistry Research

Article

(27) Woerlee, G. F.; Berends, J.; Olujic, Z.; Graauw, J. D. A comprehensive model for the pressure drop in vertical pipes and packed columns. J. Chem. Eng. 2001, 84, 367. (28) Stephan, M.; Mayinger, F. Experimental and analytical study of countercurrent flow limitation in vertical gas/liquid flows. Chem. Eng. Technol. 1992, 15, 51. (29) Nusselt, W. Die Oberflächenkondensation des Wasserdampfes. VDI Z. (1857-1968) 1916, 60, 541. (30) Ausner, I.; Hoffmann, A.; Repke, J.-U.; Wozny, G. Experimentelle und numerische Untersuchungen mehrphasiger Filmströmungen. Chem. Ing. Tech. 2005, 77, 735. (31) Paschke, S.; Repke, J.-U.; Wozny, G. Untersuchungen von Filmstromungen mittels eines neuartigen Micro particle Image Velocimetry Messverfahrens. Chem. Ing. Tech. 2008, 80, 1477. (32) Raynal, L. F.; Ben, R.; Royon-Lebeaud, A. Use of CFD for CO2 absorbers optimum design: From local scale to large industrial scale. Energy Procedia 2009, 1, 917. (33) Handbook of FLUENT 6.3, ANSYS Inc.: 2007. (34) Paschke S., Repke J-U, Wozny G. A way to analyze liquid film flows on opaque plate materials. 4th International Berlin Workshop on Transport Phenomena with Moving Boundaries, Berlin, Sept 27−28, 2007; 143.

7809

dx.doi.org/10.1021/ie500047a | Ind. Eng. Chem. Res. 2014, 53, 7797−7809