CET Volume 86 DOI: 10.3303/CET2186199 Paper Received: 6 September 2020; Revised: 8 February 2021; Accepted: 8 April 2021 Please cite this article as: Palhares P., Guirardello R., 2021, Numerical Simulation of Hold-up and Velocity Profiles for Two-phase Developed Flows in a Column, Chemical Engineering Transactions, 86, 1189-1194 DOI:10.3303/CET2186199 CHEMICAL ENGINEERING TRANSACTIONS VOL. 86,2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-84-6; ISSN 2283-9216 Numerical Simulation of Hold-Up And Velocity Profiles for Two-Phase Developed Flows in a Column Pâmela Palhares*, Reginaldo Guirardello School of Chemical Engineering, University of Campinas, Av. Albert Einstein 500, 13083-852, Campinas-SP, Brazil. pamelacoimbrap@gmail.com The purpose of this article is to simulate and analyse two-phase fluid dynamics in a bubble column, through two different numerical approaches, in a fast and effective way, to minimize project costs and increase productivity. Bubble column reactors are commonly employed in chemical processes because of their wide application area, such as simple construction and operation, high heat and mass transfer coefficients. Despite the model chosen in this article does not consider the reaction rate, bubble columns usually feature complex hydrodynamic interactions that requires a high computational effort for modelling. In this study, to simulate a bubble column reactor with continuous flow (e.g., air-water), two approaches were compared to predict axial and radial profiles for the two-phase, the finite volume and the variational method. The values of holdups and the axial velocity are in agreement with the experimental results in the literature. Key Words: fluid dynamic, hold-up, finite volume method, variational method, modelling. 1. Introduction Bubble columns are a multiphase reactors or contactors which performs interaction between gas and liquid phases with the purpose of promoting heat and mass transfer. This equipment is widely used in chemical, biochemical and petrochemical industries for having excellent heat and mass transfer coefficient, low maintenance and operating costs, simple design and operation such as the absence of moving parts. Despite its simple structure the column provides a good environment for the complex fluid dynamics interactions, which is the focus of this work. The mixing behaviour and complex flow found in bubble columns are linked to the geometry column and the operation conditions, and often can be described as the global gas holdup and liquid circulation velocity (Wu and Al-Dahhan, 2001). The gas holdup ( ) is a dimensionless parameter that expresses the volumetric gas fraction over the total volume inside the column. This parameter is of great importance in the design of bubble columns for being directly attached to the phenomena of mass transport between the gas and liquid phases. The characterization of fluid dynamics has a significant effect on the operation and performance of a bubble column. According to Hyndam et al. (1997) and Deckwer et al. (1980), the experimental results obtained by parameters investigations, strictly depend on the prevailing column regime. The flow regime is classified and maintained according to the superficial gas velocity and the column diameter. In multiphase flows, the interaction between the phases is determined by solving the momentum equation which includes interfacial forces, e.g., drag force, lift force, turbulent dispersion and virtual mass. Eulerian- Eulerian approach often uses the drag force as the predominant interfacial force to predict flow pattern and hydrodynamics properties in the bubble column. Moreover, literature has shown that although the effect of drag force is significantly higher than the others interfacial forces, it is necessary to consider the lift force to solve the equation bring more accuracy to the bubble column simulation, improving predictions of the axial liquid velocity, turbulent kinetic energy and radial gas holdup (Gupta and Roy, 2013) and (Tabib and Roy, 2008). 1189 2. Mathematic model The equations regarding the fluid dynamics in a bubble column are based on the continuity equations and momentum balance, using a Eulerian approach enabling the write for cylindrical or squared column systems with axis symmetry. The general mass and momentum equation for two-phase flow can be found in Torvik and Svendsen (1990), and Grienberger and Hofmann (1992). For this paper the fluid dynamic model considers a general heterogeneous system with two-phases in a fully developed flow with all quantities e.g., holdups ( ), radial velocity ( ), axial velocity ( ) varying only in the radial direction. The two-phase are designated here by i,j. The equations for the steady state isothermal two-phase flow in cylindrical coordinates with axisymmetric (r, z) are given by: - Mass balance: ∙ [ ∙ ∙ ∙ − ∙ ∙ ∙ ( ∙ )] = (01) For this paper the mass transfer was not considered, = 0. - Momentum balance: Momentum equation in radial and axial direction, respectively, are given by: ∙ ( ∙ ∙ ∙ ∙ ) = ∙ ∙ + ∙ + + − ∙ + ∙ ∙ ∙ ∙ ∙ 2 ∙ − ∙ ∙ 2 ∙ (02) ∙ ∙ ( ∙ ∙ ∙ ∙ ) = ∙ ∙ ∙ + − ∙ + ∙ ∙ ∙ ∙ ∙ (03) The turbulent viscosity ( ( )) was calculated using a zero order turbulence model proposed by Chen et al. (1995) and Menzel et al. (1990): = + ( ) (04) The diffusivity of one phase into another ( ) affects the distribution of the phases volume fraction ( ) inside the column. The eq. (05) shows that the diffusivity coefficient can be equal to turbulent viscosity divided by the density of a given phase (i,j). = ( ) = ( ) = ( ) = ∙ ∙ ∙ 1 + 2 ∙ 1 − (05) The momentum balances account for two interfacial forces. The Magnus force, commonly known as the transversal lift is a phenomenon in which the rotation of a body modifies the trajectory of the fluid. It is generated by velocity gradients at relative speed across the bubble diameter, forcing it to move against the liquid phase velocity gradient. = 0,5 ∙ ∙ ∙ ∙ − ∙ ∙ − + ∙ − (06) The other interfacial force used is the drag force. It is a resistance force offered by the system, slowing down the bubbles’ movement due to the flow of liquid phase near at the bubble’s boundary layer. = ∙ ∙ ∙ ( − ) (07) - Boundary conditions: = 0 = 0 (08) = 0 = 0 (09) 1190 = 0 = (10) = 0 = (11) 3. Numerical Methods 3.1 Finite Volume Method The finite volume method is a class of finite difference method where the domain is subdivided into a finite number of smaller control volumes. The main characteristics of this method is the simple and easy derivation where the equations can be interpreted in physical terms as the tendency to solve partial different equations to iterative calculations through its non-linear equations. This method was widely used by Patankar and Spalding (1972), that contributed to fill the gap in the calculation procedures for three-dimensional flow at that time. Currently is one of the most used methods to calculate fluid dynamics. In this study, Excel was chosen to run this simulation because of its simple calculation platform and easy access. 3.2 Variational Method The variational method is a numerical method that allows to obtain an accurate or estimated solution for an objective function, with unknown parameters and functions. The unknown terms are determined by the operating conditions and restrictions that are applied to the system. The principle of this method is to transform a differential equation into an integral equation, mathematically equivalent, by using the Euler- Lagrange equation (Gal-Or et al., 1972). The variational model used in this work was described in Guirardello (2019). This method is solved by GAMS software and usually is defined to maximize or minimize the objective function. The software solves the designated problem in order to deliver the optimum value. 4. Simulation To investigate the fluid dynamics inside of a column, two main components of phases distributions were predicted, axial liquid velocity and gas hold-up, one-dimensional in the radial direction. The simulation was conducted for a bubble column with 29 cm inner diameter, operating in a fully developed flow, with a co- current upward flow of liquid and gas phases (air-water). In this study no chemical reaction was considered. A zero order turbulence model was assumed to express the liquid dispersion and the fluid is considered Newtonian (Menzel et al., 1990; Chen et al., 1995). The simulation was conducted in GAMS software and Excel to calculate the turbulent flow of two-phase bubble column through variational and finite volume methods. The cases were solved using 20 intervals in the radial direction and a lift coefficient ( ) of -1. In GAMS code program was used CONOPT3 to solve nonlinear equations. The values of drag coefficient ( ) and root mean square of pressure at the wall divided by apparent mixture density were estimated to find an agreement between experimental and predicted results. 5. Results and discussion Two main local parameters were evaluated in the simulations of the fluid dynamics, the gas holdup and axial velocity. Two cases simulation were conducted in order to analyse the profiles predicted for different flow regimes by distinct superficial gas velocity. In the model used local gas holdup was calculated by using the volumetric fraction balance equation (12) and the liquid holdup (13). + = 1 (12) ∙ = ( ) ∙ (13) The axial liquid velocity was calculated by sum of the eq. (03) for gas and liquid phases, canceling the drag and lift force ( = − , = − ). It was also considered that physical properties of the liquid phase e.g., the density and viscosity, are extremely higher than the gas phase ( ≪ , ≪ ) and that radial liquid velocity is significantly lower than the axial liquid velocity ( ≪ ). − ∙ + ∙ ∙ ∙ + ∙ ∙ ∙ ∙ ∙ = 0 (14) 1191 The axial gas velocity was calculated considering the sum of momentum balance eq. (03) for the gas and liquid phases. It was also considering that radial liquid velocity is significantly lower than the axial liquid velocity, ( ≪ ). ∙ ∙ − − ∙ ∙ ∙ ∙ ∙ ∙ + ∙ − = 0 (15) The radial liquid velocity was calculated from the eq. (4) divided by , then subtracting the liquid equation from the gas equation. Also was considered that the gas phase density and viscosity are extremely lower the liquid phase. ∙ ∙ ∙ ∙ ∙ 2 ∙ − ∙ ∙ 2 ∙ − ∙ ∙ = 0.5 ∙ ∙ − ∙ ∙ − (16) For case 1, it was considered a superficial gas velocity (Ug) of 0,04 m/s and superficial liquid velocity (Ul) of 0,01 m/s, with a = 7650 cm²/s², the drag force was calculated using = 40 g.cm -3.s-1. The predictions for the main local parameters from the numerical methods used were compared with the experimental results from Yao et al. (1991) presented in Figure1. Others parameters that were calculated are presented in table 1, indicating the values distributed along the 20 intervals evidenced by the ratio distribution in the first column. The 2nd to 5th column refer to the simulation via GAMS and the 6th and 7th through Excel. The effective viscosity ( ) was equal to both methods while the radial liquid velocity had a difference less than 0,5%. The Axial and radial gas velocity had a slightly difference between the finite volume and variational method. Table 1: Comparison calculated parameters in case 1 (all units in CGS). r 0.00 21.0748 0.0000 59.3642 0.0000 59.4746 0.0000 0.73 21.1272 0.0803 58.8547 -0.2943 58.9996 -0.2945 1.45 21.2813 0.1551 57.3265 -0.5791 57.5790 -0.5719 2.18 21.5276 0.2148 55.0991 -0.8262 55.3304 -0.8112 2.90 21.8503 0.2553 52.3066 -1.0196 52.4170 -1.0005 3.63 22.2273 0.2767 49.0869 -1.1542 49.0256 -1.1344 4.35 22.6300 0.2817 45.5843 -1.2312 45.3422 -1.2134 5.08 23.0239 0.2741 41.9381 -1.2566 41.5346 -1.2425 5.80 23.3676 0.2579 38.2725 -1.2390 37.7433 -1.2292 6.53 23.6139 0.2364 34.6919 -1.1877 34.0789 -1.1819 7.25 23.7090 0.2121 31.2789 -1.1109 30.6242 -1.1084 7.98 23.5928 0.1865 28.0971 -1.0154 27.4398 -1.0151 8.70 23.1990 0.1606 25.1951 -0.9055 24.5700 -0.9063 9.43 22.4549 0.1346 22.6113 -0.7832 22.0493 -0.7841 10.15 21.2813 0.1083 20.3802 -0.6478 19.9109 -0.6477 10.88 19.5930 0.0808 18.5395 -0.4944 18.1955 -0.4917 11.60 17.2984 0.0499 17.1385 -0.3111 16.9679 -0.3026 12.33 14.2993 0.0111 16.2478 -0.0696 16.3502 -0.0471 13.05 10.4915 -0.0488 15.9597 0.3030 16.6171 0.3654 13.78 5.7644 -0.1728 16.3184 0.9673 18.6248 1.1542 14.50 0.0010 0.0000 0.0000 0.0000 0.0000 0.0000 In the case 2, it was considered a Ug of 0,08 m/s and Ul of 0,01 m/s, = 14400 cm²/s², the drag force between phases was calculated using = 40 g.cm-3.s-1 The local values for axial liquid velocity and gas holdup calculated were compared with the experimental results from Torvik and Svendsen (1990) presented in Figure2. The axial and radial velocities and the effective viscosity was calculated for these conditions are showed in Table 2. However, it only displays the results from GAMS because the values obtained are quite similar and among them the variational method was the one that showed a lightly better axial liquid velocity profile. 1192 Table 2: Calculated parameters in case 2 (all units in CGS). The values presented in Table 3 refer to the mean squared error of the main parameters calculated in the cases 1 and 2. These values are referred to experimental dada used in each case. The 2nd. and 3rd. columns refer to finite volume method, and the 3rd. and 4th. columns refer to variational method. Figure1- Results for (a) axial liquid velocity and (b) local gas holdup, comparing with experimental results from Yao et al. (1991). Figure 2- Results for (a) axial liquid velocity and (b) local gas holdup, comparing with experimental results from Torvik and Svendsen (1990). r 0.00 28.9140 76.5568 0.0000 0.0000 0.73 28.9859 75.7951 -0.3488 0.2279 1.45 29.1973 73.5100 -0.6863 0.4391 2.18 29.5353 70.1974 -0.9807 0.6071 2.90 29.9780 66.0510 -1.2132 0.7202 3.63 30.4952 61.2730 -1.3768 0.7790 4.35 31.0478 56.0747 -1.4718 0.7913 5.08 31.5881 50.6605 -1.5045 0.7682 5.80 32.0597 45.2135 -1.4849 0.7209 6.53 32.3977 39.8871 -1.4239 0.6590 7.25 32.5281 34.8032 -1.3317 0.5896 7.98 32.3687 30.0553 -1.2163 0.5170 8.70 31.8284 25.7138 -1.0831 0.4438 9.43 30.8074 21.8341 -0.9348 0.3707 10.15 29.1973 18.4652 -0.7701 0.2968 10.88 26.8811 15.6600 -0.5831 0.2194 11.60 23.7328 13.4893 -0.3588 0.1326 12.33 19.6181 12.0600 -0.0651 0.0239 13.05 14.3939 11.5390 0.3674 -0.1366 13.78 7.9083 12.1610 1.0264 -0,4195 14.50 0.0010 0.0000 0.0000 0.0000 1193 Table 3: Mean squared error obtained (all units in CGS). Case 1 0.0008 26.0817 0.0008 24.5105 Case 2 0.0021 63.0019 0.0021 52.8591 According information presented in Kantarci et al. (2005), the operation condition and column geometry presented in the 1st case is considered a homogeneous bubble regime while in the 2nd case it is a churn turbulent regime. Despite of its different flow regimes seen in Figure 1(b) and 2 (b) it has shown a good similarity on the axial liquid velocity curves between experimental dada and calculations results. Although the case 2 has a higher quadratic error. Probably it was caused by the increase in the superficial gas velocity in the column forming larger bubbles destabilizing the bubble flow due to the coalescence effect. As a result, it increases the recirculation patterns inside the column. 6. Conclusions The numerical methods proposed for a two-phase bubble column showed good results, being able to give reliable results for the two cases studied. Proving that for one-dimensional system both methods have precise results. A nonlinear problem is usually sensitive to initial estimates, therefore variational method program used the values of liquid hold-ups from finite volume method to run calculations. That may be the reason why some of the calculated values are similar in both methods. For future studies some improvement on the programing code will be necessary. This study also proves that Excel, can be used to calculate fluid dynamics problems without an expensive optimization software with a proper mathematical method, minimizing project cost. For this study all physical properties were considered constant, calculation should be revised in future problems for chemical reaction or if the physical properties change with pressure. Acknowledgments The authors gratefully acknowledge the financial support of CNPq - Conselho Nacional de Desenvolvimento Científico e Tecnológico, Brazil. Reference Chen Z., Zheng C., Feng. Y., 1995, Modelling of three-phase fluidized beds based on local bubble characteristics measurements, Chemical Engineering Science, 50, 231–236. Deckwer W.D., Louisi Y., Zaidi A., Ralek M., 1980, Hydrodynamic properties of Fisher-Tropsch slurry process, Ind. Eng. Chem. Process Des. Dev., 19, 699-708. Gal-Or B, Wheihs D., 1972, Variational analisys of high mass transfer rates from spherical particles boundary layers, Int. J. Mass Heat Transfer, 15,2027-2044. Grienberger J., Hofmann H., 1992, Investigations and modelling of bubble columns, Chemical Engineering Science, 47, 2215–2220. Guirardello R., 2019, A Numerical Resolution of Fluid Dynamic Problems Using a Saddle Point Variation Formulation, Chemical Engineering Transactions, 74, 1015-1020. Gupta A., Roy S., 2013, Euler–Euler simulation of bubbly flow in a rectangular bubble column: experimental validation with Radioactive Particle Tracking, Chemical Engineering Journal, 225, 818-836 Kantarci N., Boral F., Ulgen K.O., 2005, Bubble column reactor, Process Biochemistry, 40, 2263-2283. Hyndam C.L., Larachi F.,Guy C.,1997, Understanding gas-phase hydrodynamics in bubble column: a convective model based on kinetic theory, Chemical engineering science, 52, 63-77. Patankar, S. V & Spaldíng, D.B., 1972, A Calculation Procedure for Heat, Mass and Momentum Transfer in Three-Dimenssional Parabolic Flow, Int. J. Heat Mass Transfer, 15, 1787. Pourtousi M., Sahu J.N., Ganesan P.,2014 effect of interfacial forces and turbulence models on predicting flow pattern inside the bubble column, Chemical engineering and processing, 75,38-47. Tabib M.V., Roy S.A., Joshi J.B., 2008, CFD simulation of bubble column—an analysis of interphase forces and turbulence models, Chemical Engineering Journal, 139, 589-614. Torvik R., Svendsen H.F., 1990, Modelling of slurry reactors. A fundamental approach, Chemical Engineering Science, 45, 2325–2332. Yao B.P., Zheng C., Gasche H. E., Hofmann H.,1991, Bubble behaviour and flow structure of bubble columns, Chemical Engineer Process, 29, 65-75. Wu Y., Al-Dahhan M.H., 2001 Prediction of axial liquid velocity profile in bubble column, Chemical engineering Science, 56, 1127-1130. 1194