部分论文openFoam资料Three-dimensional mathematical modeling of dispersed two-phase flow

Computers and Chemical Engineering 32 (2008) 32243237 Contents lists available at ScienceDirect Computers and Chemical Engineering journal homepage: Three-dimensional mathematical modeling of dispersed two-phase fl ow using class method of population balance in bubble columns Rachid Bannaria, Fouzi Kerdoussb, Brahim Selmaa, Abdelfettah Bannaria, Pierre Proulxa, aDepartment of Chemical Engineering, Universit de Sherbrooke, Sherbrooke, QC J1K 2R1, Canada bEnvironnement Canada, Place Vincent Massey, 351, Boulevard St-Joseph, Gatineau, QC K1A 0H3, Canada a r t i c l ei n f o Article history: Received 26 October 2007 Received in revised form 6 May 2008 Accepted 28 May 2008 Available online 8 June 2008 Keywords: Population balance equation Coalescence and break-up Interfacial area Methods of classes Bubble size Computational fl uid dynamics (CFD) a b s t r a c t Computationalfl uiddynamics(CFD)simulationsofbubblecolumnshavereceivedrecentlymuchattention and several multiphase models have been developed, tested, and validated through comparison with experimental data. In this work, we propose a model for two-phase fl ows at high phase fractions. The inter-phase forces (drag, lift and virtual mass) with different closure terms are used and coupled with a classes method (CM) for population balance. This in order to predict bubbles size distribution in the column which results of break-up and coalescence of bubbles. Since these mechanisms result greatly of turbulence, a dispersed k ? turbulent model is used. Theresultsarecomparedtoexperimentaldataavailablefromtheliteratureusingameanbubblediame- ter approach and CM approach and the appropriate formulations for inter-phase forces in order to predict the fl ow are highlighted. The above models are implemented using the open source package OpenFoam. Crown Copyright 2008 Published by Elsevier Ltd. All rights reserved. 1. Introduction Bubble columns are today widely used in many industrial pro- cesses like chemical, pharmaceutical and petrochemical because of their capability of achieving high heat and mass transfer rates with lowenergyinput.Thecomputationalfl uiddynamics(CFD)methods usingmodernhighspeedcomputingcapabilitiesisausefultoolthat is now often used in the industry for scaling-up or design of many types of reactors, but a lot of effort is to be made in order to develop models that will take into account all the complexity of bubble columnstypereactors.Soundengineeringjudgementanddetailed understanding of the underlying physics of the phenomena must be used in order to develop models that can adequately describe the observed behavior of these reactors. In bubble columns, gas phase exists as a dispersed bubble phase in a continuous liquid phase. Up to the recent years, the use of CFD modeling for bubble columns has been very limited compared to that in single-phase or dilute multiphase systems. Most work published had focused on developing closure models for interfacial bubbleliquid forces andonestimatingbubble-inducedturbulence.Again,untilrecently, relatively little attention has been given to the bubble size distri- bution problem although it is an important design parameter and Corresponding author. E-mail address: ulxusherbrooke.ca (P. Proulx). can infl uence signifi cantly the result of gasliquid mass transfer equipment. Bubble size distribution depends extensively on col- umn geometry, operating conditions, physical properties of phases and sparger type, and the last few years have seen more and more researchers addressing this complex problem. There are essentially two approaches for the numerical calcula- tion of multiphase fl ow, namely the EulerLagrange method which considers the bubbles as individual entities tracked using trajec- tory equations (Webb, Que, Lapin Pfl eger, Gomes, Gilbert, Buwa Kumar Lehr Venneker, Derksen, WangandWang,2007).Chen,Dudukovic, and Sanyal (2005a,b) concluded in their comparison with different bubblebreakupsandcoalescenceclosuresthatthechoiceofbubble breakup and coalescence closure does not have a signifi cant infl u- enceontheresultsaslongasthemagnitudeofbreakupisincreased tenfold. In the present work, the model of Luo and Svendsen (1996) is used for its simplicity, a predefi ned daughter bubble size distri- bution is not needed to predict breakup rate for bubbles ad thus avoids some of the problems of the original model. The daughter size distribution can be directly calculated from the model. In the present work all the inter-phase forces (drag, lift and virtual masse) (Mudde Buwa Bel FDhila Buwa, Deo, Sanyal, Marchisio,Fox,Bhole,Joshi, Rusche, 2002). Transport equations for the phase fraction as well as the turbulent kinetic energy and its dissipation rate are solved only once per time step/iteration at the end of the sequence. 2. Previous work Modeling interfacial forces remains an open question in the numerical simulation of bubbly fl ow. It represents interaction forces between the dispersed and continuous phase. When motion is non-uniform, bubbles are accelerated relatively to the liquid. The concept of drag is extended to include various non-drag forces such as the so-called virtual mass force and lateral lift force. An adequate description of this closure law must be used to predict accurately the behavior of bubbly fl ow. Sustained efforts have been made to understand the complex fl ow fi elds in a bubble column. Among early work, several authors neglected the effect of virtual mass (Deb Roy, Majumdar, Joshi, 1981, 1983; Kumar, Devanathan, Moslemian, Mudde Chen et al., 2005a,b; Kerdouss, Bannari, Bel FDhila Lance Lathouwers,1999)thattheliftforcehasastronginfl uenceonphase fraction distribution and it is widely acknowledged that the stan- dard lift model which contains an empirical lift coeffi cient gives widely different results in different types of fl ow. In recent years, there has been strenuous efforts to reformulate and understand interfacial drag force, virtual mass and lift force (Ranade, 1992; Sokolichin Becker, Sokolichin, Phanikumar Sokolichin Krishna, van Baten, Buwa Behzadi, Issa Buwa et al. 2006; Zhang, Deen, Mudde Pfl eger Troshko ?i(v) = xi+1 v xi+1 xi (41) Death in class i due to coalescence DiC= 36fi2 g ?2d3 i n ? k=0 a(xi,xk) fk d3 k (42) Birth in class i due to breakup BiB= 6g ? n ? k=i m(xk)b(xk)?i,k fk d3 k (43) Death in class i due to breakup DiB= 6gfi ?d3 i b(xi)(44) and ?i,k= ? xi xi1 v xi1 xi xi1 p(v,xk)dv + ? xi+1 xi xi+1 v xi+1 xi p(v,xk)dv(45) Using Gaussian quadrature integration, the above integrals are approximated with ?i,k? 5 ? j=1 (1 + Wj)3 (j + 1)2P2 5(Wj) p(xi xi1 2 (1 + Wj) xi1,xk) + 5 ? j=1 (1 + Wj)2(1 Wj) (j + 1)2P2 5(Wj) p(xi+1 xi 2 (1 + Wj) xi,xk) (46) Pnis a Legendre polynomial which can be constructed using the three term recurrence relations: Pn= (2n 1)xPn1 (n 1)Pn2 n ,P0= 1, P1= x(47) Wjis the weighting function related to the orthogonal polyno- mials, see Table 4. 3230R. Bannari et al. / Computers and Chemical Engineering 32 (2008) 32243237 Table 4 Values of weighting function used in Gaussian quadrature integration W1 ? 35+270 63 W2 ? 35270 63 W30.0 W4= W2?352 70 63 W5= W1?35+2 70 63 3.7. Bubble breakup Bubble break-up is described as a consequence of a collision of a bubble with a turbulent eddy (Prince Kerdouss et al., 2008): B(vj,vi) = 18k1lgfj 11?b8/11d3 j ? ? d2 j ?1/3 ?(8/11,tm) ?(8/11,b) + 2b3/11(?(5/11,tm)?(5/11,b) +b6/11(?(2/11,tm) ?(2/11,b) (49) where b = 12cf? 1?c?2/3d5/3 j ;tm= b(?min/dj)11/3(50) 1? 2.05 and k1? 0.924. At high Reynolds numbers, terms with tmare taken equals to zero as tm? (Sanyal et al., 2005). By this defi nition of the break-up kernel, terms b(v?) and p(v,v?) can be written as b(v?) = 1 m(v?) ? v? 0 B(v?,v)dv(51) This integral is calculated numerically in the discretized popu- lation equation. p(v,v?) = B(v?,v) ?v? 0 B(v?,v)dv (52) 3.8. Bubble breakup and coalescence Coalescence rates a(vi,vj) are usually written as the product of the collision rate ?ij and coalescence effi ciency Pc(Hagesather, Jakobsen, Hjarbo, ?ij= di/dj;uij= (u2 i + u2 j) 1/2; ui= 1/2(?dweusei)1/3(56) 4. Boundary conditions All walls are treated as non-slip boundaries with standard wall function (Versteeg and Malalasekera, 1995, p. 200). The gas fl ow rate at the sparger is defi ned via inlet velocity type boundary con- dition with the gas volume fraction equal to unity. The bubble size at the gas inlet depends on the sparger design which is beyond the scope of this work and an uniform bubble size at the inlet of 5 mm was taken here based on the study of Buwa et al. (2006). At the liq- uid surface, a gas zone is added at the free surface of water in order to prevent liquid escape from the column, only gas is allowed to escape (g= 1, l= 0). The boundary condition for ? and k at inlet are as follows: k = 3 2(UinletI) 2 (57) where I is the turbulence intensity ? = C3/4 ? k3/2 0.07L (58) with L the characteristic length of the equipment or pipe radius. 5. Numerical solution In this work, the Open source Field Operation and Manipu- lation (Open FOAM) C+ libraries are used. They are based on a tensorial approach developed fi rst by Weller, Tabor, Jasak, and Fureby (1998). It is freely available and open source, licensed under the GNU General Public Licence. OpenFOAM is supplied with numerous pre-confi gured solvers, utilities and libraries, It is open, not only in terms of source code, but also in its struc- ture and hierarchical design, so that its solvers, utilities and libraries are fully extensible. It uses fi nite volume method to solve systems of partial differential equations described on any 3D unstructured mesh of polyhedral cells. New solvers and utili- ties can be created by the user with some pre-requisite knowledge of the underlying method, physics and programming techniques involved. The aim in OpenFOAM is to offer an unrestricted choice to user andcompletefreedomtochoosefromawideselectionofinterpola- tion schemes. In this work, linear Gaussian integration for gradient R. Bannari et al. / Computers and Chemical Engineering 32 (2008) 322432373231 Table 5 Two-phase numerical solution procedure 1. Calculate the inter-phase force (closure model A or B) 2. Solve the momentum Eq. (4) 3. PISO-Loop 4. Correct the substantive derivatives 5. Solve k ?Eqs. (18) and (19) 6. Solve the gEq. (16) 7. Solve the PBE equations (if needed) operators is chosen. It is based on summing value on cell faces, which must be interpolated from cell centers. For divergence terms (ex. convection term is the momentum equation) the Gauss lin- ear scheme for discretization is used, it calculates interpolations based on fl uxes. For the fi rst time derivative, the choice is the Euler implicit fi rst order (Hrvoje, 1996; Rusche, 2002; Open Foam user guide 2006, 2007). The solution of the resulting equations is made with the segregated technique, where within an iterative cycle, each set of algebraic equations is solved in turn using an appropriate solver until convergence is obtained. A special treat- ment is required in order to establish necessary inter-equation coupling (Hrvoje, 1996; Rusche, 2002). Solution of the pressure equationprovidescorrectionsforupdatingpressure,fl uxandveloc- ities, so that continuity is satisfi ed. Since two phases continuity equations are used, a mixture continuity form (Eq. (6) (Hill, 1998; Rusche, 2002). PISO (Pressure Implicit with Splitting of Operators) algorithm proposed by Issa is used to handle pressurevelocity coupling. It is based on procedures described by Weller (2002), where a pressure equation based on volumetric continuity equa- tionissolvedtocorrectamomentumpredictedbecauseofpressure change. This is in a correction loop which consists of an implicit momentum predictor followed by a series of pressure solutions and explicit velocities corrections. Iterative methods are used to solve the system of algebraic equations created by discretization. Diagonal dominance which guarantees convergence, is taken into account (Jasak, 1996) and methods choice is dependent of matrix structure. Incomplete Cholesky preconditioned Conjugate Gradi- ent (ICCG) solver is used for symmetric matrices. The method is described in detail by Jacobs (1980). The solver for asymmet- ric matrices is Bi-CGSTAB by Van Der Vorst (1992). For pressure, ICCG with tolerance equal to 1e 8 and BICCG for velocities, dis- sipation, turbulent kinetic energy and phase fraction is used. The sequence of operation for the solution procedure is summarized in Table 5. The calculations are carried out for a rectangular bubble col- umn similar to that used by Pfl eger et al. (1999), and used by Buwa and Ranade (2002, 2006). It is 0.2 m width 1.2 m height 0.05m depth. The water level is 0.45m (H/W = 2.25) and 0.9m (H/W = 4.5). Superfi cial gas velocities used are 0.14 and 0.73cm/s. Pseudo structured grids are used, 92% structured with 6100 cells for H/W ratio of 2.25 and 93.5% structured with 30,232 cells for H/W ratio of 4.5. Buwa and Ranade (2002) carried out calculations usingdifferentconfi gurationsoflocallyaeratedspargersandfound that results are not sensitive to the different sparger representa- tions. Therefore, sparger through which gas was introduced into the column is modeled as area covered by sparging holes (18mm 6mm in the present case). Typical grid distribution and sparger representation are shown in Fig. 1. A time step of 1.0 103s was used. For the PBE, the bubble size is divided into n = 2r + 1classes, with n odd in order to have symmetry. As described before a dis- tribution on pivotal grid points xiwith xi+1= sxiand s 1 is used, whereirefertotheclassiwithi n.Withtheassumptionofspher- ical bubbles, we can write (4/3)?(di+1/2)3= (4s/3)?(di/2)3. s is calculated to ensure that dn= d2r+1= dmaxand dr= dmean Fig. 1. Typical grids used in present work with H/W ratio of 4.5 and 2.25. This gives the following relation: di= s(ir1/3) dmeanands = ? dmax dmean ?3/r . The bubble mean diameter is set to 5 mm and maximal diam- eter to 10mm as in Buwa and Ranade (2002) and Buwa et al. (2006).Table 6 gives the values of s used with the number of classes used with CM. At the inlet the mean diameter is used. In this work, the same velocity Ug is used for all size groups. This simplifi cation is used to reduce the gas phase momentum equations to a single equation with a velocity Ug. Further work is underway to take into account the difference between velocities of different classes but is not included here. 6. Results and discussion ExperimentsbyBuwaandRanade(2002)andBuwaetal.(2006) show that when gas is introduced into a column fi lled with liquid, Table 6 Values of s and the different number of classes used with CM Number of classes7111525 Value of r35712 Value of s21.51571.34591.1892 3232R. Bannari et al. / Computers and Chemical Engineering 32 (2008) 32243237 Fig. 2. Experimental and modeling (gas volume fraction, 0 black and 0.05 white) snapshots of meandering bubble plume at a superfi cial gas velocity of 0.73cm/s (H/W = 2.25) by Buwa et al. (2006). gas bubbles formed at sparger holes rise upwards in pool of liquid exhibiting different length and time scales. Experimental instanta- neoussnapshotsofoscillatingbubbleplume,computationalresults published by Buwa et al. (2006) and the predictions of the cur- rent model are shown in Figs. 2 and 3. Predicted instantaneous gas volume fraction distribution show that meandering motion of the bubbleplumeiscapturedinsatisfactoryqualitativeagreementwith experiment. For a quantitative comparison of the simulated results Fig. 3. Gas volume fraction (0 blue and 0.05 red), and liquid velocity snapshots of meandering bubble plume at a superfi cial gas velocity of 0.73cm/s (H/W =2.25) using the current model. (For interpretation of the references to colour in this fi gure legend, the reader is referred to the web version of the article.) Fig.4. Voidagefl uctuationtimeseriesobtainedfromthecurrentmodel(withsuper- fi cial gas velocity: 0.14cm/s; H/W: 2.25). Fig. 5. Time-averaged vertical liquid velocity profi les predicted using closure A and B with one mean bubble diameter and the measurement by Pfl eger et al. (1999) at an H/W ratio of 2.25 and a superfi cial gas velocity of 0.14cm/s (Y = 0.37m). with experimental measurements, time-averaged fl ow properties (vertical liquid velocity and gas hold-up) are compared with mea- surements using two types of closure model, namely model A and model B in Table 2 with different classes of bubble (010 mm) or using mean bubble diameter (5 mm), as used by Buwa et al. (2002, 2006). Voidage fl uctuation time series obtained using the cur- rent model, based on the two phases fl uid methodology of Weller (2002),atasuperfi cialgasvelocityof0.14cm/s,anH/W ratioof2.25 and recorded at X = 0.1m, Y = 0.25m and Z = 0.025m, is shown in Fig. 4. It agrees well with experimentally recorded voidage fl uc- Fig. 6. Time-averaged gas hold-up predicted using closure models A and B with one mean diameter and the measurement and predited results by Buwa et al. (2006) at an H/W ratio of 2.25 and a superfi cial gas velocity of 0.14cm/s (Y = 0.37 m). R. Bannari et al. / Computers and Chemical Engineering 32 (2008) 322432373233 Fig.7. Time-averagedverticalliquidvelocityprofi lespredictedusingcurrentmodel with closure A (mean diameter), the CM (between 7 and 25 classes) and the mea- surementbyPfl egeretal.(1999)atanH/W ratioof2.25andasuperfi cialgasvelocity of 0.14cm/s (Y = 0.37m). Fig.8. Time-averagedverticalliquidvelocityprofi lespredictedusingcurrentmodel with closure B (mean diameter), the CM (between 7 and 25 classes) and the mea- surementbyPfl egeretal.(1999)atanH/W ratioof2.25andasuperfi cialgasvelocity of 0.14cm/s (Y = 0.37m). tuation time series using conductivity probe published by Buwa et al. (2006). Indeed, the number of low-frequency oscillations corre- spondingtothemotiono