部分论文openFoam资料KlostermannJ Numerical simulation of a single rising bubble by VOF with surface compression

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN FLUIDS Int. J. Numer. Meth. Fluids (2012) Published online in Wiley Online Library ( DOI: 10.1002/fl d.3692 Numerical simulation of a single rising bubble by VOF with surface compression J. Klostermann*, K. Schaake and R. Schwarze Institute of Mechanics and Fluid Dynamics, TU Bergakademie Freiberg, Freiberg, Germany SUMMARY The capability of the direct volume of fl uid method for describing the surface dynamics of a free two- dimensional rising bubble is evaluated using quantities of a recently published benchmark. The model equations are implemented in the open source computational fl uid dynamics library OpenFOAM. Here, a main ingredient of the numerical method is the so-called surface compression that corrects the fl uxes near the interface between two phases. The application of this method with respect to two test cases of a bench- mark is considered in the main part. The test cases differ in physical properties, thus in different surface tension effects. The quantities centre of mass position, circularity and rise velocity are tracked over time and compared with the ones given in the benchmark. For test case one, where surface tension effects are more pre-eminent, deviations from the benchmark results become more obvious. However, the fl ow features are still within reasonable range. Nevertheless, for test case two, which has higher density and viscosity ratios and above all a lower infl uence of the surface tension force, good agreement compared with the benchmark reference results is achieved. This paper demonstrates the good capabilities of the direct volume of fl uid method with surface compression with regard to the preservation of sharp interfaces, boundedness, mass conservation and low computational time. Some limitation regarding the occurrence of parasitic currents, bad pressure jump prediction and bad grid convergence have been observed. With these restrictions in mind, the method is suitable for the simulation of similar two-phase fl ow confi gurations. Copyright 2012 John Wiley Revised 16 April 2012; Accepted 30 April 2012 KEY WORDS:VOF; surface tension force; parasitic/spurious currents; direct method; free surface; rising bubble benchmark 1. INTRODUCTION Interfacial fl ows of two or more immiscible fl uids are widely used in many industrial processes, for example, in casting operations, food processing, fi bre coating, emulsifi cation, spraying and so on. The kinematics and dynamics of the fl uid interfaces play an important role in many of these processes. For example, the interface between the covering slag and the liquid steel in a continuous casting mould can strongly affect the quality of the fi nal product, when interfacial instabilities are induced by the fl ow 1. Fast, robust and accurate numerical models are necessary for the investigation of the interfacial fl uid dynamics by means of computational fl uid dynamics (CFD) on fi xed computing grids. Promi- nent approaches are marker and cell 2, volume of fl uid method (VOF) 3 and level set method (LS) 4,5. A detailed discussion of these approaches is beyond the scope of this paper; interested readers are referred to review papers, for example, by Rider and Kothe 6, Scardovelli and Zaleski 7 or Sethian and Smereka 8. *Correspondence to: J. Klostermann, Institute of Mechanics and Fluid Dynamics, TU Bergakademie Freiberg, Lampadiusstr. 4, 09599 Freiberg, Germany. E-mail: jens.klostermannimfd.tu-freiberg.de Copyright 2012 John Wiley for the simulations later on, the CrankNicolson time scheme is used. Copyright 2012 John Wiley here, the divergence term is rewritten with Gausss theorem: ? D ?r ?.E n/ D ? X f .b ni/f(28) 2.2.3. Solution algorithm. The overall solution cycle is given in Algorithm 1. Here, the distribu- tion of the volume fraction fi eld ( ?-subcyle at step 4) is calculated ahead of the well-known PISO algorithm 26 (step 8) that updates the velocity and pressure fi elds. The key element of the ?-subcyle is the solution of Equation (24). To guarantee a sharp and bounded solution, a numerical scheme designed for the multi-dimensional advection equation should be employed (e.g. 27,28). For this study, the explicit multidimensional universal limiter with explicit solution (MULES) implemented in OpenFOAM20, 21 is used to integrate the ?-equation Equation (24) explicitly. Therefore, Equations (23) and (24) are solved together with the mass fl ux ?f?i ?f?iWD ?fE uf? E S D .?1/f? C.1?/?2?f? D ?.?1?2/fC.?2/f?(29) for a number of N?corrcorrector steps. During this correction, ?rstays constant, see Algorithm 1. Finally, the mass fl ux is calculated according to Equation (24) as ?f? D N?sub X iD1 ti ?t ?f?i(30) with tiD ?t=N?sub. Equations (4), (6), (25) and (30) are iterated inside the so-called ?-sub cycle. This sub cycle consists of solving the equations N?sub-times (sub time-stepping; here, also ?r is updated, see in Algorithm 1) and then averaging the mass fl ux. During the PISO loop, the mass fl ux ?f? in the advection term is maintained constant. Once the pressurevelocity equations are solved, the volumetric fl ux ? is updated. The overall numerical procedure is sketched in Algorithm 1. 2.2.4. Discretisation schemes and used parameters. The applied discretisation schemes and the parameters of the numerical model are summarised in Table I. For convenience, the corresponding terminology of OpenFOAMis also given. The transient term in the momentum equation, Equation (14), is discretised by the CrankNicolson scheme 29. Different interpolation schemes are applied to determine cell face values in the convective terms. For the convective part in the momentum equation, a total variation diminishing scheme with the fl ux limiter function .r/ D max0,min.2?r,1/? is used. The smoothness parameter r is defi ned by the ratio of consecutive gra- dients. The van Leer 30 scheme with the fl ux limiter function .r/ D .rCjrj/=.1Cjrj/ is used for thediscretisationofther?.E u?/terminEquation(13).Here,thefl uxisboundedbetween0and1.For the second convective part r ? .E ur.1?/?/ in Equation (13), the so-called interfaceCompression scheme is used. Here, the limiter function reads .?P, ?N/ D min ? max 1?max ?q 1?4?P?.1?p/, p 1?4?N?.1?N/ ? ,0 ,1 ? (31) which is bounded between 0 and 1, with the fl ux ?P at the evaluated cell and the fl ux ?Nat the neighbour cell. In our study, we included the optional momentum predictor step that brought small but noticeable improvements in the results. This will be discussed later. Copyright 2012 John Wiley CDS, central differencing scheme. ?The symbol ? stands for an arbitrary scalar or vector. 3. NUMERICAL SET-UP OF THE SIMULATIONS 3.1. Computational domain On the basis of the benchmark defi nition of Hysing et al. 9, a two-dimensional computational domain with an aspect ratio x W y D 1 W 2 is employed, see Figure 1. The bubble is initially centred at .x,y/ D .0.5,0.5/ with rb0 D 0.25 as the initial radius. The viscosity and density of fl uid 2 (bubble ?2 ) are smaller than those of the surrounding fl uid 1 (?1). The domain is fully enclosed by no-slip walls (E u D .0,0/) at the top and the bottom and free slip walls (uxD 0 and r? f E u D 0) on the left and the right. The gravity vector E g points towards the bottom of the domain. Hysing et al. 9 distinguished two different set-ups of the numerical experiment: ellipsoidal bub- ble (TC1) and skirted bubble (TC2). They differ by the density ratio ?1=?2and the ratio of the viscosity ?1=?2 between both phases (index 1 for heavy and index 2 for light fl uid). Physical parameters of TC1 and TC2 are given in Table II. All parameters in Table II are made dimensionless with characteristic scales for the length L D 2rb0, for the time t D L=Ugand for the rising velocity UgD pg2r b0. The Reynolds number Re, the Etvs number Eo and the capillary number Ca are defi ned as Re D ?1UgL ?1 ,Eo D ?1U2 gL ? ,Ca D ?1Ug ? D Eo Re (32) Thus, the surface tension force is more prominent for TC1 in comparison with TC2, see Table II. The simulations of TC1 and TC2 are performed with four different grid spacings h D 1=40, 80, 160, 320?, respectively. The simulation time is t fi nalD 3 with a grid size-dependent time step of ?t D h=2. Copyright 2012 John Wiley in practice, the region is spread over a few cells so that the surface position can vary within this region. To pinpoint the surface, we chose a value of ? D 0.5. 3.2.3. Maximum rise velocity Vmaxwith corresponding incidence time. The maximum rise velocity is deduced from the rise velocity V , which is evaluated for each time step as follows. V D R ?1?2 ? v dA R ?1?2 ? dA (35) Here, v is the velocity in y-direction in an individual cell in the computational mesh. (a) domain after patching(b) domain after relaxation Figure 2. Distribution of both phases in the computational domain of the coarsest grid (h D 1=40), yellow corresponds to ? D 1, blue to ? D 0 and green to 0 ? 1, situation after patching and relaxation. (a) Domain after patching and (b) relaxation. Copyright 2012 John Wiley no cells with 0 ? 0. In detail, the development of maxjupj shows small dependence on the grid spacing and on c?. In the case of ?r 0 (c?D 1), a convergent behaviour for a resolution fi ner than h 6 1=40 is noticeable. Note that the parasitic current is pro- portional to the kinetic energy in the system. To obey energy conservation, the kinetic energy and thus the parasitic currents should approach zero because of the inclusion of viscous forces. A snapshot of the pressure fi eld for TC1 under zero gravity condition at t D 3 is depicted in Figure 5. The pressure is normalised by the pressure jump across the surface given by the Young Laplace equation ?p D 2 ?=rb0resulting in p? D p=?p D p rb0=.2?/ . In the case of a static bubble, the normalised pressure should read in ?1: p? D 0 and ?2: p? D 1 with a sharp pressure jump at the surface. This is not the case because the values range from p? D 0.83 to p? D 0.70 (in ?2) for the mesh resolution ranging from h D 1=40 to h D 1=320, respectively. Furthermore, at the interface, some unphysical spikes (overshots and undershots) in the pressure fi eld can be found. To reduce these spikes, we varied the interpolation scheme (least squares and cell-limited linear schemes) 20,21 for the r? term, which gave no improvements. The resulting errors of the zero gravity simulations for TC1 and TC2 on the coarsest and the fi nest grids are given in Table IV. The order of magnitude of Capremains nearly constant in both test cases. No grid convergence was found. Therefore, Cap? 10?4is determined as a measure to estimate the infl uence of the parasitic currents in our simulations. However, the parasitic currents only have a marginal infl uence on the centre of mass position because this error is much smaller than the nominal displacement length of the parasitic currents, (a)(b) Figure 4. Infl uence of compression fl ux ?ron the parasitic currents. (a) c?D 0 and (b) c?D 1. Copyright 2012 John Wiley thus, we use only TP2D to refer both of them in this test case. Obviously, the VOF solutions on the coarsest and the fi nest grids are very close to each other, but there is a small but noticeable difference to the fi nal bubble shape of the reference results, see Figure 7 (left). A detailed view at a typical section reveals that congruence between the VOF simulations and the reference results seems to be unreachable, not even with much fi ner grid resolutions. Figures 810 illustrate the temporal development of the quantities centre of mass position Yc, circularity C and rising velocity V from which the benchmark quantities are deduced. Again, the results of the VOF simulations for the four-grid resolutions h D 1=40, 80, 160, 320? are compared with the corresponding results of the codes TP2D and FreeLIFE given by Hysing et al. 9. The temporal evolution of Ycin our simulations is nearly identical to the reference results, see Figure 8 (left), but the close-up showing the range t D 2.5?3 in Figure 8 (right) shows that complete convergence towards the benchmark results is again not achieved. On the contrary, the temporal evo- lution of C in Figure 9 (left) and V in Figure 10 (left) shows more apparent differences between our fi ndings and the corresponding benchmark data. Here, the close-ups near the minimum circularity Cminin Figure 9 (right) and the maximum rising velocity Vmaxin Figure 10 (right) reveal that the |The reference results have the following grid resolutions h D 1=640 (TP2D) and h D 1=160 (FreeLIFE) respectively. Copyright 2012 John Wiley hence, all simulations are carried out with the additional momentum predictor step. 4.3. TC2 Figure 15 shows typical results for the rising bubble in TC2 at three different stages (t D 0, t D 1.5 and t D 3). The distortion of the interface in TC2 is considerably larger than in TC1. The bubble is deformed to a skirted bubble and eventually forms thin fi laments with smaller satellite bubbles that are detached. Compared with TC1, lower surface tension forces indicated by the higher Eo and Ca numbers as well as higher viscosity and density ratios are the reasons for the formation of a skirted bubble. The correct simulation of these physics is especially challenging as there is Copyright 2012 John Wiley even so, convergence, that is a grid-independent solution, is again not achieved. Contrary to the fi ndings for TC1, the evolution of C in Figure 18 (left) and V in Figure 19 (left) of our simulations are also closer to the corresponding results of FreeLIFE. Please note the differences in the benchmark results: in TP2D, there is a separation of satellite bubbles, whereas in FreeLIFE, the bubble remains skirted with elongated fi laments. This difference between TP2D and FreeLIFE is documented by the very different temporal evolution of C. In our simulations, the rising bubble also remains skirted, although the separation of the satellite bubbles is found close to the fi nal situation at t D 3. The close-ups near the minimum circularity Cminin Figure 18 (right) and the maximum rising velocities Vmax1and Vmax2in Figure 19 (right) confi rm that the results of our simulations approximate the results of FreeLIFE, although complete convergence is again not achieved. Copyright 2012 John Wiley 162167. 2. Harlow F, Welch J. Numerical calculation of time-dependent viscous incompressible fl ow of fl uid with free surface. Physics of Fluids 1965; 8(12):21822189. 3. Hirt CW, Nichols BD. Volume of fl uid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics 1981; 39:201225. 4. Osher S, Sethian JA. Fronts propagating with curvature-dependent speedalgorithms based on HamiltonJacobi formulations. Journal of Computational Physics 1988; 79(1):1249. 5. Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase fl ow. Journal of Computational Physics 1994; 114(1):146159. Copyright 2012 John Wiley 141(2):112152. 7. Scardovelli R, Zaleski S. Direct numerical simulation of free-surface and interfacial fl ow. Annual Review of Fluid Mechanics Jan 1999; 31(1):567603. 8. Sethian J, Smereka P. Level set methods for fl uid interfaces. Annual Review of Fluid Mechanics 2003; 35:341372. 9. Hysing S, Turek S, Kuzmin D, Parolini N, Burman E, Ganesan S, Tobiska L. Quantitative benchmark computations of two-dimensional bubble dynamics. International Journal for Numerical Methods in Fluids 2008; 60:12591288. 10. Ferziger JH, Peri c M. Computational Methods for Fluid Dynamics. Springer: Berlin, Germany, 2002. 11. Yeoh GH, Barber T. Assessment of interface capturing methods in computational fl uid dynamics (CFD) codesa case study. The Journal of Computational Multiphase Flows Jun 2009; 1(2):201215. 12. trubelj L, Tiselj I, Mavko B. Simulations of free surface fl ows with implementation of surface tension and interface sharpening in the two-fl uid model. International Journal of Heat and Fluid Flow 2009; 30:741750. 13. Ubbink O. Numerical prediction of two fl uid systems with sharp interfaces. PhD Thesis, Imperial College, University of London, 1997. 14. Muzaferija S, Peric M. Computation of free-surface fl ows using the fi nite-volume method and moving grids. Numerical Heat Transfer Part B-Fundamentals December 1997; 32(4):369384. 15. YoungsD.Time-dependentMulti-material fl owwithLargeFluidDistortion.AcademicPressInc: (London)Ltd,1982. 227(11):53615384. 19. zkan F, Wrner M, Wenka A, Soyhan H. Critical evaluation of CFD codes for interfacial simulation of bubble- train fl ow in a narrow channel. International Journal for Numerical Methods in Fluids 2007; 55(6):537564. DOI: 10.1002/fl d.1468. 21. OpenCFD Limited. User guide openfoam 1.6. Technical Report. http:/www.opencfd.co.uk/openfoam;2009. 22. Brackbill JU, Kothe DB, Zemach C. A continuum method for modeling surface tension. Journal of Computational Physics 1992; 100(2):335354. 23. Rusche H. Computational fl uid dynamics of dispersed two-phase fl ows at high phase fractions. PhD Thesis, Imperial College of Science, Technology and Medicine, London, England, 2002. 24. de M, P B R. Study and numerical simulation of sediment transport in free-surface fl ow. PhD Thesis