Direct numerical simulation of rising bubble interaction ... - Springer Link
Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148 www.springerlink.com/content/1738-494x
DOI 10.1007/s12206-012-0819-3
Direct numerical simulation of rising bubble interaction with free surface using level contour reconstruction method† Seungwon Shin* School of Mechanical and System Design Engineering, Hongik University, Seoul, 151-791, Korea (Manuscript Received March 1, 2012; Revised May 6, 2012; Accepted May 23, 2012) ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Abstract Bubble rising phenomenon is widely found in many engineering applications including stream generators in power plants. Many experimental and numerical researches have been performed to predict the dynamic behavior of the bubble rising process. Most simulations so far have focused on evolution behavior of the rising bubble itself. Rising bubbles could penetrate through the top free surface which makes the problem much more complicated in addition to the curvature effect on the interface. With top surface, satellite droplets may be generated near breakup region and entrained to rising velocity field. In this paper, the effect of top free surface on rising bubble has been numerically investigated. High-order tetra-marching level contour reconstruction method (LCRM), which is a hybridization of fronttracking and level-set methods has been used to track the gas-liquid interface explicitly. Effects of free surface on the interfacial shape and behavior after breakup are studied with different size and depth of the initial bubble. It has been found that the initial diameter of the bubble is the dominant factor controlling the satellite droplet formation. Keywords: Numerical simulation; Multiphase flow; Bubble rising; Top free surface; Front tracking; Level set ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1. Introduction Many industrial heat generation power plant systems utilize phase changing multiphase flow, considering its large heat capacity in regard to energy transportation. Recently, the combined heat and power (CHP) plant has drawn quite an attention for the interest of generating both electrical power and heat at the same time. It can increase the overall system performance considerably. Mid or small sized CHP plants have been constructed for relatively small business buildings, and sometimes it requires steam generation for the sanitization purpose. Thus, the heat recovery steam generator (HRSG) has been installed in CHP for the necessary steam demand. Steam drum at the top of the HRSG would be one of the key components to capture and separate the steam from the boiling water. To increase the steam generation rate, it is very important to understand the process inside the steam drum correctly. Despite its importance, experimental analysis of rising bubbles inside the drum has many difficulties due to small length and time scale associated with the bubble rising process. Numerical simulation can be very efficient for identifying the underlying physical behavior. *
Corresponding author. Tel.: +82 320 3038, Fax.: +82 320 7003 E-mail address: [email protected] Recommended by Associate Editor Gihun Son © KSME & Springer 2012 †
The bubble rising problem has been studied extensively since it is related to key features for various engineering applications, including bubble columns and vapor generation devices, to name a few [1, 2]. The bubble rising phenomenon itself also provides a very fundamental aspect in regard to benchmarking numerical techniques for multiphase flows [3, 4]. Lately, numerical methods using fixed Eulerian grid with additional advection scheme to track the interface motion, e.g. VOF [3], level set [4], and front tracking [2], have been widely used. Most numerical simulations so far have concentrated on single bubble rise and focused on investigating the shape and terminal velocity of the rising bubble. In case of the rising bubble’s confronting a top free surface, interface evolution will be quite different from the single bubble rising process. The rising bubble will penetrate through the top free surface and satellite droplets might be generated above the free interface. It is important to describe the behavior of the bubble and the top free surface correctly after collision, since the formation of the satellite droplet on top of the free surface will have considerable influence on the efficiency of the steam separation process. Since there will be complicated interface evolution related to the bubble rise with free surface, an accurate and reliable numerical method is required to track the interface motion precisely and model complex topological changes appropriately at the same time. Nowadays, hybridization of existing
3142
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
(a)
Fig. 1. Basic concept of the level contour reconstruction method.
numerical techniques such as VOF, level set, and front tracking to overcome inherent problems associated with each method have become popular. Level set methods [5, 6] based on finite element method and unstructured grid have also been proposed. Some examples are CLSVOF (level set+VOF) [79] originally developed by Sussman and Puckett, particle level set (level set+Lagrangian marker particle) from Enright et al. [10], and Aulisa et al. [11] who combine markers and VOF method. Level contour reconstruction method (LCRM) [1214] which combines the advantage of front tracking and level set methods can be also considered as hybrid type methods. The level contour reconstruction method is basically a front tracking type method which explicitly tracks implicitly connected individual interface elements. Since the Lagrangian interface can be reconstructed using Eulerian function field (e.g., distance function) which can be computed directly from given interface location, merging and pinch off of interfaces can be dealt with naturally and automatically as in the level set method. In this study, the rise of two-dimensional axisymmetric vapor bubble with free surface will be simulated using the level contour reconstruction method. Focus has been made to identify the controlling factors for the bubble dynamics with the help of hybrid numerical scheme. Effects of free surface on the interfacial shape and behavior after breakup are studied with different size and depth of the initial bubble. Effect from the ambient pressure on the bubble behavior has also been identified and a parametric study has been performed to find dominant variables on steam drum design.
2. Numerical formulation 2.1 Interface tracking method Level contour reconstruction method [12-14] is based on the fact that the interface can be represented by both Lagrangian interface and certain contour of Eulerian field variables. As can be seen from Fig. 1, indicator function, I, which is basically Heaviside type function or distance function, φ, can be generated with existing Lagrangian interfacial information. At the same time, we can reconstruct Lagrangian interface
(b) Fig. 2. Tetra marching procedure for the two-dimensional level contour reconstruction method: (a) directional ambiguity in drawing contour line; (b) triangular cell orientation for tetra marching process.
elements by drawing contour lines with value 0.5 of indicator function field or 0 of distance function field. The motion of the evolving interface will be tracked explicitly with implicitly connected Lagrangian elements, and the interface elements will be periodically reconstructed using the Eulerian distance function field. During the original reconstruction procedure [12], the interface has been reconstructed using linear interpolation from the given distance function field. We found that linear interpolation generates a continuous but not a smooth interface after reconstruction. High order LCRM [13] has been introduced to increase the accuracy and smoothness of the reconstructed interface. In general, the reconstruction procedure can be performed with same rectangular Eulerian grid structure for flow computation. With this rectangular shaped grid, there can be more than two surfaces inside the single reconstruction cell. Special formulation of interface reconstruction is required for rectangular shaped reconstruction domain to avoid the ambiguity in connecting surface elements as shown in Fig. 2(a). Interface reconstruction of the LCRM doesn't have to coincide with a rectangular grid cell. In case of using a triangular cell for interface reconstruction, there can be a single possibility for drawing contour line around edges as in Fig. 2(b). We try to use given Cartesian rectangular cell for fluid computation as our basic reconstruction cell and subdivide this into two triangular cells. Tetrahedral reconstruction cell has been reoriented as shown in Fig. 2(b) to guarantee connectivity of reconstructed interface elements. More detailed description of tetramarching procedure can be found in Ref. [15]. For bubble rise with top surface, severe topological changes can be anticipated. With conventional contour reconstruction, interface elements will be merged together if the two separate elements become close enough within reconstruction grid resolution. With high-order tetra-marching reconstruction,
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
3143
symmetric version. Detailed explanation of the hybrid surface tension formulation can be found in Ref. [14]. The local surface tension force at the interface, F, can be described by the following equation: F = σκ H ∇I .
(a)
(b)
Fig. 3. Topological control of the LCRM: (a) controlling of the merging distance; (b) distributing weight function for topological control of the iterface.
interface elements can be closer than reconstruction cell size. With front tracking capability associated with the LCRM, the proximity of the interface elements can be easily identified as depicted in Fig. 3(a). If two interface elements are closer than a user defined constant (dl in Fig. 3(a)), then the distance function value of the reconstruction cell near the interface element is weighted with a certain value to force the merging or breakup of the interfaces (Fig. 3(b)). In this way we can control the subgrid distance for merging and breakup as in the conventional front tracking method. More detailed numerical structure of the LCRM can be found in Refs. [12-15]. 2.2 Numerical formulations For the current problem, we have used an axisymmetric version of the LCRM to improve efficiency of computations, since it is still difficult to conduct full simulation in threedimensional domain. Single field formulation of the governing equations in non-conservative axisymmetric form for multifluid motion can be expressed as [16]: 1 ∂ru ∂v + =0 r ∂r ∂z
∂u ∂u ∂p ∂ ∂u ∂u +u +v =− + 2µ ∂r ∂z ∂r ∂r ∂r ∂t
(1)
ρ
+2µ
∂ u ∂ ∂v ∂u + µ + + F ⋅ ir ∂r r ∂z ∂r ∂z
∂v ∂v ∂p 1 ∂ ∂v ∂u ∂v +u +v =− + + µr t r z ∂ ∂ ∂ ∂z r ∂r ∂r ∂z , ∂ ∂v + 2 µ + F ⋅ iz − ρ g z ∂z ∂z
(4)
Here, σ is the surface tension coefficient (assumed constant here), I is the indicator function, a Heaviside function which varies from zero to one near the interface. Indicator function also has been utilized to describe material property at each phase. The Heaviside function can be calculated using distance function, φ, which can be computed directly from given interfacial information. Detailed procedure for computing distance function can also be found in Ref. [14]. The hybrid form of the curvature field, κH, on the Eulerian grid can be obtained using
κH =
FL ⋅ G
σG ⋅G
.
(5)
In this equation the force, FL , is the front tracking formulation for the surface tension force which has been modified to account for axisymmetric component as: FL =
∫ σκ f
fnfδ
( x − x f )ds − nrr
(6)
where nf is the unit normal to the interface, xf = x(s,t) is a parameterization of the interface, κf is the local mean curvature calculated on the interface, and δ(x-xf) is a twodimensional Dirac distribution that is non-zero only when x = xf. The first term on the right hand side of Eq. (6) represents the regular surface tension force term in plane 2D simulation. The second term shows the added curvature effect from circumferential rotation of axisymmetric component. In the second term, nr is interface normal and r is the interface position in r direction, respectively. The distribution, G, is geometric information computed directly on the interface and then distributed onto an Eulerian grid. G=
∫ n δ ( x − x )ds f
(2)
f
f
(7)
The interface is advected in a Lagrangian fashion by integrating
ρ
(3)
Here, u and v are the velocity in r and z directions, p is the pressure, gz is the gravitational acceleration in z direction, and F is the interfacial surface tension force term. We have used the hybrid formulation [14] with some modification for axi-
dx f dt
=V
(8)
where V is the interface velocity vector interpolated at xf. More detailed solution procedure and discretization of the governing equations can be found in Refs. [12-14]. Chorin’s projection method [17] has been used to solve the fluid velocity and pressure. Staggered mesh (MAC method) of
3144
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
Fig. 5. Benchmarking test for comparison of the rising velocity with grid convergence.
Fig. 4. Schematic description of the simulation domain for axisymmetric geometry.
Harlow and Welch [18] has been used for spatial discretization with the forward Euler method of time integration. The pressure is located at the cell centers, while the r and z components of velocity are located at the respective faces. All spatial derivatives except the convective term are approximated by standard second-order centered differences. The convective term is discretized by using a second-order ENO procedure [19]. The detailed solution procedure and discretization of the governing equations can be found in Refs. [1214]. Fig. 6. Mass conservation of the single bubble rise test with different grid resolutions.
3. Results and discussion 3.1 Benchmarking test of single bubble rising simulation To verify the accuracy of the axisymmetric code and find optimum grid resolution, we have considered a twodimensional bubble rise problem. Single bubble rise in a liquid pool has been studied widely with various numerical techniques to verify its accuracy [1, 2, 7, 20, 21]. Spherical bubble released in a liquid pool will reach steady rising velocity after some developing period. The rising velocity of the bubble is determined by initial size of the bubble and property ratio between vapor and ambient liquid. The computational domain for axisymmetric bubble rise simulation is shown in Fig. 4. For this benchmarking test, we have not included the top free surface. Initially, a spherical droplet with unit length (D = 1) is placed at the lower bottom of the simulation domain. The left side of the simulation domain is the axis of symmetry. Sufficiently large domain size has been used for axial and radial length of the simulation domain to neglect the effect from the boundary. Box size of 3D × 12D has been used for the simulation. Open condition has been given to upper boundary, and bottom boundary was treated as a stationary wall. The governing Eqs. (1), (2), and
(3) have been non-dimensionalized using the scale of length as initial bubble diameter, D, velocity as Uo = (gD)1/2, and pressure as ρLUo2 (measured from ambient system pressure Patm). The problem can then be characterized by the Reynolds, Re, and Weber, We, numbers as well as the property ratios which can be defined as follows: Re =
ρLU o D ρ U 2D , We = L o µL σ
ρ∗ =
ρG ∗ µG ,µ = . µL ρL
(9) (10)
Here, ρ and µ represent density and viscosity, respectively. Subscript L is for the ambient liquid properties and G for the gas bubble. For single bubble rise test, Re = 8.66 and We = 1.26 have been used to compare simulation results from Ryskin and Leal [1], which were computed with body fitted grid. Both density and viscosity ratio have been set to 0.001. Fig. 5 shows non-dimensional rising velocity with increasing grid resolution. With increasing number of grids, we can see that
3145
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
Table 1. Material properties at various ambient pressure for numerical simulation. Patm [atm] 1 3
the simulation result is converging to the value of Ryskin and Leal [1]. In Fig. 6, mass conservation of the rising bubble with different grid resolution has also been plotted. Relatively large amount of mass has been lost with lowest resolution. This mass loss will be the major reason for the deviation of the rising velocity at the lowest resolution in Fig. 5. Interface plot in the early stage of the rising process is shown in Fig. 7 with three different grid resolutions at simultaneous time. We can clearly see that the interface evolution is converging with increasing grid resolution. Considering the simulation results from Figs. 5-7, a grid resolution of 120 × 720 has been chosen as the working resolution for current numerical simulation tests. 3.2 Bubble rising simulation with top free surface The bubbles generated from the evaporation column will be collected in the steam drum and then separated from free surface. The separated steam should be exhausted without any droplets for the efficient operation of the HRSG system. The dominant forces determining the physical behavior of the interface motion would be the gravity, surface tension, and inertia. We have chosen initial diameter (D in Fig. 4), rising distance (H in Fig. 4), and ambient pressure as the major controlling parameters. The material properties of the bubble (vapor in Fig. 4) and water (liquid in Fig. 4) will be a function of the ambient pressure. Considering the general operational condition of the HRSG system, we have selected a few possible values of the ambient pressure for current parametric study. We have summarized the simulation condition for the various input pressure in Table 1. We have simulated a two-dimensional axisymmetric bubble rise problem with top free surface as described in Fig. 4.
5
10
20
ρL [kg/m ]
958.4
942.7
914.9
886.6
849.1
ρG [kg/m3]
0.5969
1.142
2.699
5.206
10.17
µL×104 [N-s/m]
2.819
2.308
1.795
1.498
1.257
µG×10 [N-s/m]
0.1227
0.1298
0.1407
0.1504
0.1617
σ [N/m]
0.05891
0.05484
0.04824
0.04208
0.03467
4
Fig. 7. Interface profile of the single bubble rise simulation for three difference grid resolutions.
2
Simulation condition is exactly the same as the single bubble rise case in section 3.1 except top free surface presence. The radial length of the simulation domain stays the same as single bubble rising simulation, but the height of the simulation domain changes according the rising distance, H. The properties of steam water at 5 atm, typical operating condition of a small HRSG system, has been chosen as the reference case. Property ratio becomes ρ* = 0.00295 for the density and µ* = 0.07839 for the viscosity, respectively. Initially, spherical bubble has been placed at (0, 3.5D) position. The box size of 3D × 11D has been used for the reference case. The resulting interface evolution has been plotted in Fig. 8. The rising bubble becomes elongated in longitudinal direction to form an ellipsoidal bubble, then hits the top free surface. The thin filamentary interface between bubble and top surface breaks up to form a satellite droplet at the top of the free surface. For the bubble rise with top surface, the focus has been to identify the satellite droplet formation after bubble merging to the top free surface. The formation and behavior of the satellite droplet after impact would be the critical factor for proper operation of the steam drum in the HRSG system. The best possible scenario would be no droplet generation since an entrained droplet at the upper steam collector would damage the system performance. The droplet might explode with oversaturation temperature or condense the surrounding steam to form larger droplet, which ultimately hinders steam flow. Even if there is a droplet formation after impact, sufficiently small size of the droplet could be entrained due to the magnitude of upward steam velocity from phase change process at the free surface. So the size of the satellite droplet would be the important output parameter determining the performance of steam drum system. In addition to the size of the satellite droplet, its velocity would also determine the destination of the droplet. With large velocity, the detached droplet would travel further to reach the steam collecting area. To analyze the behavior of the bubble with top free surface more relevantly for the HRSG system, we chose our output parameters as detaching time, time required for the rising bubble to interact with top free surface, diameter of detached droplet, maximum velocity of the detached droplet. We tested the grid convergence of the bubble rise simulation with top free surface. The non-dimensional grid size for the reference case (see Fig. 8) is 0.025, which has been identi-
3146
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
Fig. 8. Interface evolution of the single bubble rise with top surface at several time steps.
(a)
(b) Fig. 9. Grid convergence test for single bubble rise with top surface: (a) convergence of detaching time vs. grid size; (b) convergence of diameter of detached drop vs. grid size.
fied as the working resolution for the single bubble rise test. We further decreased or increased the grid size to investigate its effect on simulation results. The non-dimensional detaching time of the satellite droplet and non-dimensional diameter of the detached droplet vs. grid size is plotted in Fig. 9. As can be seen there, the working grid size of 0.025 is again appropriate for the bubble rise with top surface. The merging distance, dl, in Fig. 3 has been set to 0.025, single grid size, for
the all simulation considered. As discussed before, active controlling parameters for the droplet formation after merging of the rising bubble with top surface are rising distance, initial diameter of the bubble, and ambient pressure for the material properties. To identify the individual effect from each input variable of the rising distance, initial diameter of the bubble, and ambient pressure, we fixed the other two variables and varied only one variable of interest from the reference case. In Fig. 10, the effect of rising distance (H), ambient pressure (Patm), and initial bubble diameter (D) on non-dimensional diameter of detached droplet is shown. As can be seen, rising distance and ambient pressure have less effect than initial droplet on the size of detached droplet. With increasing distance, H, and ambient pressure, Patm, the size of the detached droplet is increasing a little but relatively small amount. The initial size of the bubble has the dominant effect on the size of the satellite droplet. Initial bubble, which is less than 1 mm or greater than 8 mm in diameter, doesn’t form any satellite droplet. It can be readily expected that there will be no detached droplet with small initial bubble. The inertial force from the small rising bubble is not enough to make a detached surface after impact. In Fig. 13, we have shown interface evolution with large initial diameter of the bubble (D = 30 mm). The rising bubble tends to elongate further with drag force in horizontal direction with increasing initial volume. So the lower portion of the bubble becomes hollow to make a spherical cap shaped interface. The bottom surface will ultimately hit the topside of the bubble to make a donut shaped interface as in Fig. 13. This ring type interface will rise further to merge with top free surface but does not have enough inertia to from a satellite droplet. It is interesting to note that there is a maximum of the detached droplet size in the middle of the considered size distribution of the initial bubble. The detaching time is exactly proportional to rising distance, H, as expected as described in Fig. 11. The bubble rising velocity is converging to a constant very quickly as can be seen from Fig. 5. The non-dimensional detaching time has very weak effect from the ambient pressure, similar to the size of the detached droplet. Again, initial diameter of the bubble is the major controlling parameter for the detaching time. The detaching time is increasing with larger initial bubble size. The maximum velocity of the detached droplet is depicted in Fig. 12 with each controlling parameter. Again, rising distance and ambient pressure have an insignificant effect on the maximum velocity of the detached droplet. With changing ambient pressure, there is a minimal value near the reference pressure of 5 atm even though relative changes are fairly low. Initial diameter of the bubble has a major effect on the final velocity of the detached droplet. In this case, there is a minimum value of the maximum velocity as opposed to the result of detached droplet size. It can be readily inferred from Fig. 10 that the size of the detached droplet is inversely proportional to the maximum velocity.
3147
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
(a)
(b)
(c)
Fig. 10. Effects of input variables on the non-dimensional diameter of detached drop: (a) effect of distance to the top surface; (b) effect of ambient pressure; (c) effect of initial bubble diameter.
(a)
(b)
(c)
Fig. 11. Effects of input variables on the non-dimensional detaching time: (a) effect of distance to the top surface; (b) effect of ambient pressure; (c) effect of initial bubble diameter.
(a)
(b)
(c)
Fig. 12. Effects of input variables on the non-dimensional maximum velocity of detached drop: (a) effect of distance to the top surface; (b) effect of ambient pressure; (c) effect of initial bubble diameter.
4. Conclusions The bubble rise problem with top free surface has been analyzed numerically with level contour reconstruction method. The LCRM is a hybrid type interface method for multiphase flows which combines some essential features of front tracking and level set techniques. By adopting high order tetramarching reconstruction, we can maintain excellent interface fidelity and perfectly interconnected interfacial elements with
complex topological changes. A benchmarking test of single bubble rising has been performed to show the accuracy of axisymmetric version of the LCRM. The results shows good convergence with the existing body fitted numerical solution. For the bubble rise with top free surface, the focus has been to identify the formation of the satellite droplet. The effect from the rising distance, ambient pressure, and initial diameter of the bubble on the size and velocity of the detached droplet and detaching time has been identified. The influence of the initial size of the bubble is dominant for the interface evolution with
3148
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
the top free surface. The size of the detached droplet has maximum value on specific initial bubble diameter as opposed to the maximum velocity of a detached droplet showing minimum value. The interface generally undergoes severe transition of the topological evolution. More quantitative research of multiple bubbles interaction with top free surface, especially in 3D, is necessary and currently underway.
Acknowledgment This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (20090071479, 20100016889).
References [1] G. Ryskin and L. G. Leal, Numerical simulation of freeboundary problems in fluid mechanics. Part 2, J. Fluid Mech., 148 (1984) 19-35. [2] A. Esmaeeli and G. Tryggvason, Direct numerical simulation of bubbly flows. Part 1. Low Reynolds number arrays, J. Fluid Mech., 377 (1998) 313-345. [3] R. Scardovelli and S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Ann. Rev. Fluid Mech., 31 (1999) 567-603. [4] S. Osher and R. P. Fedkiw, Level set methods: An overview and some recent results, J. Comput. Phys., 169 (2001) 463502. [5] H. G. Choi, A numerical study on the entry problem of an elastic body using a level set method and a monolithic formulation, Journal of Mechanical Science and Technology (2012) (in press). [6] H. G. Choi, A least square weighted residual method for level set formulation, Int. J. Numerical Methods in Fluids, 68 (2012) 887-904. [7] M. Sussman and E. G. Puckett, A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows, J. Comput. Phys., 162 (2000) 301-337. [8] G. Son and N. Hur, A coupled level set method and volume of fluid method for the buoyancy-driven motion of fluid particles, Int. J. of Numerical Heat Transfer, 42 (2002) 523-542. [9] J. Choi and G. Son, Numerical study of droplet motion in a microchannel with different contact angles, Journal of Mechanical Science and Technology, 22 (2008) 2590-2599. [10] D. Enright, R. Fedkiw, J. Ferziger and I. Mitchell, A hybrid particle level set method for improved interface capturing, J. Comput. Phys., 183 (2002) 83-116.
[11] E. Aulisa, S. Manservisi and R. Scardovelli, A mixed markers and volume-of-fluid method for the reconstruction and advection of interfaces in two-phase and free-boundary flows, J. Comput. Phys., 188 (2003) 611-639. [12] S. Shin and D. Juric, Modeling three-dimensional multiphase flow using a level contour reconstruction method for front tracking without connectivity, J. Comput. Phys., 180 (2002) 427-470. [13] S. Shin and D. Juric, High order level contour reconstruction method, Journal of Mechanical Science and Technology, 21 (2007) 311-326. [14] S. Shin and D. Juric, A hybrid interface method for threedimensional multiphase flows based on front tracking and level set technique, Int. J. Numerical Methods in Fluids, 60 (7) (2009) 753-778. [15] I. Yoon and S. Shin, Tetra-marching procedure for high order level contour reconstruction method, Advances in fluid mechanics VIII, WIT press (2010) 507-518. [16] J. Han and G. Tryggvason, Secondary breakup of axisymmetric liquid drops. I. Acceleration by a constant body force, Physics of Fluids, 11 (12) (1999) 3650-3667. [17] J. Chorin, Numerical solution of the Navier-Stokes equations, Math. Comput. 22 (1968) 745-762. [18] F. H. Harlow and J. E. Welch, Numerical calculation of time dependent viscous incompressible flow of fluid with free surface, Phys. Fluids. 8 (1965) 2182-2189. [19] C. W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes II, J. Comp. Phys. 83 (1989) 32-78. [20] M. Sussman, K. M. Smith, M. Y. Hussaini, M. Ohta and R. Zhi-Wei, A sharp interface method for incompressible twophase flows, J. Comp. Phys. 221 (2007) 469-505. [21] J. Hua, J. F. Stene and P. Lin, Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method, J. Comp. Phys. 227 (2008) 3358-3382.
Seungwon Shin received his B.S. and M.S. degrees in Mechanical Engineering from Seoul National University, Korea, in 1995 and 1998, respectively. He then received his Ph.D. from Georgia Tech in 2002. Dr. Shin is currently a Professor at the School of Mechanical and System Design Engineering at Hongik University in Seoul, Korea. His research interests include computational fluid dynamics, multiphase flow, surface tension effect, phase change process.
DOI 10.1007/s12206-012-0819-3
Direct numerical simulation of rising bubble interaction with free surface using level contour reconstruction method† Seungwon Shin* School of Mechanical and System Design Engineering, Hongik University, Seoul, 151-791, Korea (Manuscript Received March 1, 2012; Revised May 6, 2012; Accepted May 23, 2012) ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Abstract Bubble rising phenomenon is widely found in many engineering applications including stream generators in power plants. Many experimental and numerical researches have been performed to predict the dynamic behavior of the bubble rising process. Most simulations so far have focused on evolution behavior of the rising bubble itself. Rising bubbles could penetrate through the top free surface which makes the problem much more complicated in addition to the curvature effect on the interface. With top surface, satellite droplets may be generated near breakup region and entrained to rising velocity field. In this paper, the effect of top free surface on rising bubble has been numerically investigated. High-order tetra-marching level contour reconstruction method (LCRM), which is a hybridization of fronttracking and level-set methods has been used to track the gas-liquid interface explicitly. Effects of free surface on the interfacial shape and behavior after breakup are studied with different size and depth of the initial bubble. It has been found that the initial diameter of the bubble is the dominant factor controlling the satellite droplet formation. Keywords: Numerical simulation; Multiphase flow; Bubble rising; Top free surface; Front tracking; Level set ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1. Introduction Many industrial heat generation power plant systems utilize phase changing multiphase flow, considering its large heat capacity in regard to energy transportation. Recently, the combined heat and power (CHP) plant has drawn quite an attention for the interest of generating both electrical power and heat at the same time. It can increase the overall system performance considerably. Mid or small sized CHP plants have been constructed for relatively small business buildings, and sometimes it requires steam generation for the sanitization purpose. Thus, the heat recovery steam generator (HRSG) has been installed in CHP for the necessary steam demand. Steam drum at the top of the HRSG would be one of the key components to capture and separate the steam from the boiling water. To increase the steam generation rate, it is very important to understand the process inside the steam drum correctly. Despite its importance, experimental analysis of rising bubbles inside the drum has many difficulties due to small length and time scale associated with the bubble rising process. Numerical simulation can be very efficient for identifying the underlying physical behavior. *
Corresponding author. Tel.: +82 320 3038, Fax.: +82 320 7003 E-mail address: [email protected] Recommended by Associate Editor Gihun Son © KSME & Springer 2012 †
The bubble rising problem has been studied extensively since it is related to key features for various engineering applications, including bubble columns and vapor generation devices, to name a few [1, 2]. The bubble rising phenomenon itself also provides a very fundamental aspect in regard to benchmarking numerical techniques for multiphase flows [3, 4]. Lately, numerical methods using fixed Eulerian grid with additional advection scheme to track the interface motion, e.g. VOF [3], level set [4], and front tracking [2], have been widely used. Most numerical simulations so far have concentrated on single bubble rise and focused on investigating the shape and terminal velocity of the rising bubble. In case of the rising bubble’s confronting a top free surface, interface evolution will be quite different from the single bubble rising process. The rising bubble will penetrate through the top free surface and satellite droplets might be generated above the free interface. It is important to describe the behavior of the bubble and the top free surface correctly after collision, since the formation of the satellite droplet on top of the free surface will have considerable influence on the efficiency of the steam separation process. Since there will be complicated interface evolution related to the bubble rise with free surface, an accurate and reliable numerical method is required to track the interface motion precisely and model complex topological changes appropriately at the same time. Nowadays, hybridization of existing
3142
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
(a)
Fig. 1. Basic concept of the level contour reconstruction method.
numerical techniques such as VOF, level set, and front tracking to overcome inherent problems associated with each method have become popular. Level set methods [5, 6] based on finite element method and unstructured grid have also been proposed. Some examples are CLSVOF (level set+VOF) [79] originally developed by Sussman and Puckett, particle level set (level set+Lagrangian marker particle) from Enright et al. [10], and Aulisa et al. [11] who combine markers and VOF method. Level contour reconstruction method (LCRM) [1214] which combines the advantage of front tracking and level set methods can be also considered as hybrid type methods. The level contour reconstruction method is basically a front tracking type method which explicitly tracks implicitly connected individual interface elements. Since the Lagrangian interface can be reconstructed using Eulerian function field (e.g., distance function) which can be computed directly from given interface location, merging and pinch off of interfaces can be dealt with naturally and automatically as in the level set method. In this study, the rise of two-dimensional axisymmetric vapor bubble with free surface will be simulated using the level contour reconstruction method. Focus has been made to identify the controlling factors for the bubble dynamics with the help of hybrid numerical scheme. Effects of free surface on the interfacial shape and behavior after breakup are studied with different size and depth of the initial bubble. Effect from the ambient pressure on the bubble behavior has also been identified and a parametric study has been performed to find dominant variables on steam drum design.
2. Numerical formulation 2.1 Interface tracking method Level contour reconstruction method [12-14] is based on the fact that the interface can be represented by both Lagrangian interface and certain contour of Eulerian field variables. As can be seen from Fig. 1, indicator function, I, which is basically Heaviside type function or distance function, φ, can be generated with existing Lagrangian interfacial information. At the same time, we can reconstruct Lagrangian interface
(b) Fig. 2. Tetra marching procedure for the two-dimensional level contour reconstruction method: (a) directional ambiguity in drawing contour line; (b) triangular cell orientation for tetra marching process.
elements by drawing contour lines with value 0.5 of indicator function field or 0 of distance function field. The motion of the evolving interface will be tracked explicitly with implicitly connected Lagrangian elements, and the interface elements will be periodically reconstructed using the Eulerian distance function field. During the original reconstruction procedure [12], the interface has been reconstructed using linear interpolation from the given distance function field. We found that linear interpolation generates a continuous but not a smooth interface after reconstruction. High order LCRM [13] has been introduced to increase the accuracy and smoothness of the reconstructed interface. In general, the reconstruction procedure can be performed with same rectangular Eulerian grid structure for flow computation. With this rectangular shaped grid, there can be more than two surfaces inside the single reconstruction cell. Special formulation of interface reconstruction is required for rectangular shaped reconstruction domain to avoid the ambiguity in connecting surface elements as shown in Fig. 2(a). Interface reconstruction of the LCRM doesn't have to coincide with a rectangular grid cell. In case of using a triangular cell for interface reconstruction, there can be a single possibility for drawing contour line around edges as in Fig. 2(b). We try to use given Cartesian rectangular cell for fluid computation as our basic reconstruction cell and subdivide this into two triangular cells. Tetrahedral reconstruction cell has been reoriented as shown in Fig. 2(b) to guarantee connectivity of reconstructed interface elements. More detailed description of tetramarching procedure can be found in Ref. [15]. For bubble rise with top surface, severe topological changes can be anticipated. With conventional contour reconstruction, interface elements will be merged together if the two separate elements become close enough within reconstruction grid resolution. With high-order tetra-marching reconstruction,
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
3143
symmetric version. Detailed explanation of the hybrid surface tension formulation can be found in Ref. [14]. The local surface tension force at the interface, F, can be described by the following equation: F = σκ H ∇I .
(a)
(b)
Fig. 3. Topological control of the LCRM: (a) controlling of the merging distance; (b) distributing weight function for topological control of the iterface.
interface elements can be closer than reconstruction cell size. With front tracking capability associated with the LCRM, the proximity of the interface elements can be easily identified as depicted in Fig. 3(a). If two interface elements are closer than a user defined constant (dl in Fig. 3(a)), then the distance function value of the reconstruction cell near the interface element is weighted with a certain value to force the merging or breakup of the interfaces (Fig. 3(b)). In this way we can control the subgrid distance for merging and breakup as in the conventional front tracking method. More detailed numerical structure of the LCRM can be found in Refs. [12-15]. 2.2 Numerical formulations For the current problem, we have used an axisymmetric version of the LCRM to improve efficiency of computations, since it is still difficult to conduct full simulation in threedimensional domain. Single field formulation of the governing equations in non-conservative axisymmetric form for multifluid motion can be expressed as [16]: 1 ∂ru ∂v + =0 r ∂r ∂z
∂u ∂u ∂p ∂ ∂u ∂u +u +v =− + 2µ ∂r ∂z ∂r ∂r ∂r ∂t
(1)
ρ
+2µ
∂ u ∂ ∂v ∂u + µ + + F ⋅ ir ∂r r ∂z ∂r ∂z
∂v ∂v ∂p 1 ∂ ∂v ∂u ∂v +u +v =− + + µr t r z ∂ ∂ ∂ ∂z r ∂r ∂r ∂z , ∂ ∂v + 2 µ + F ⋅ iz − ρ g z ∂z ∂z
(4)
Here, σ is the surface tension coefficient (assumed constant here), I is the indicator function, a Heaviside function which varies from zero to one near the interface. Indicator function also has been utilized to describe material property at each phase. The Heaviside function can be calculated using distance function, φ, which can be computed directly from given interfacial information. Detailed procedure for computing distance function can also be found in Ref. [14]. The hybrid form of the curvature field, κH, on the Eulerian grid can be obtained using
κH =
FL ⋅ G
σG ⋅G
.
(5)
In this equation the force, FL , is the front tracking formulation for the surface tension force which has been modified to account for axisymmetric component as: FL =
∫ σκ f
fnfδ
( x − x f )ds − nrr
(6)
where nf is the unit normal to the interface, xf = x(s,t) is a parameterization of the interface, κf is the local mean curvature calculated on the interface, and δ(x-xf) is a twodimensional Dirac distribution that is non-zero only when x = xf. The first term on the right hand side of Eq. (6) represents the regular surface tension force term in plane 2D simulation. The second term shows the added curvature effect from circumferential rotation of axisymmetric component. In the second term, nr is interface normal and r is the interface position in r direction, respectively. The distribution, G, is geometric information computed directly on the interface and then distributed onto an Eulerian grid. G=
∫ n δ ( x − x )ds f
(2)
f
f
(7)
The interface is advected in a Lagrangian fashion by integrating
ρ
(3)
Here, u and v are the velocity in r and z directions, p is the pressure, gz is the gravitational acceleration in z direction, and F is the interfacial surface tension force term. We have used the hybrid formulation [14] with some modification for axi-
dx f dt
=V
(8)
where V is the interface velocity vector interpolated at xf. More detailed solution procedure and discretization of the governing equations can be found in Refs. [12-14]. Chorin’s projection method [17] has been used to solve the fluid velocity and pressure. Staggered mesh (MAC method) of
3144
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
Fig. 5. Benchmarking test for comparison of the rising velocity with grid convergence.
Fig. 4. Schematic description of the simulation domain for axisymmetric geometry.
Harlow and Welch [18] has been used for spatial discretization with the forward Euler method of time integration. The pressure is located at the cell centers, while the r and z components of velocity are located at the respective faces. All spatial derivatives except the convective term are approximated by standard second-order centered differences. The convective term is discretized by using a second-order ENO procedure [19]. The detailed solution procedure and discretization of the governing equations can be found in Refs. [1214]. Fig. 6. Mass conservation of the single bubble rise test with different grid resolutions.
3. Results and discussion 3.1 Benchmarking test of single bubble rising simulation To verify the accuracy of the axisymmetric code and find optimum grid resolution, we have considered a twodimensional bubble rise problem. Single bubble rise in a liquid pool has been studied widely with various numerical techniques to verify its accuracy [1, 2, 7, 20, 21]. Spherical bubble released in a liquid pool will reach steady rising velocity after some developing period. The rising velocity of the bubble is determined by initial size of the bubble and property ratio between vapor and ambient liquid. The computational domain for axisymmetric bubble rise simulation is shown in Fig. 4. For this benchmarking test, we have not included the top free surface. Initially, a spherical droplet with unit length (D = 1) is placed at the lower bottom of the simulation domain. The left side of the simulation domain is the axis of symmetry. Sufficiently large domain size has been used for axial and radial length of the simulation domain to neglect the effect from the boundary. Box size of 3D × 12D has been used for the simulation. Open condition has been given to upper boundary, and bottom boundary was treated as a stationary wall. The governing Eqs. (1), (2), and
(3) have been non-dimensionalized using the scale of length as initial bubble diameter, D, velocity as Uo = (gD)1/2, and pressure as ρLUo2 (measured from ambient system pressure Patm). The problem can then be characterized by the Reynolds, Re, and Weber, We, numbers as well as the property ratios which can be defined as follows: Re =
ρLU o D ρ U 2D , We = L o µL σ
ρ∗ =
ρG ∗ µG ,µ = . µL ρL
(9) (10)
Here, ρ and µ represent density and viscosity, respectively. Subscript L is for the ambient liquid properties and G for the gas bubble. For single bubble rise test, Re = 8.66 and We = 1.26 have been used to compare simulation results from Ryskin and Leal [1], which were computed with body fitted grid. Both density and viscosity ratio have been set to 0.001. Fig. 5 shows non-dimensional rising velocity with increasing grid resolution. With increasing number of grids, we can see that
3145
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
Table 1. Material properties at various ambient pressure for numerical simulation. Patm [atm] 1 3
the simulation result is converging to the value of Ryskin and Leal [1]. In Fig. 6, mass conservation of the rising bubble with different grid resolution has also been plotted. Relatively large amount of mass has been lost with lowest resolution. This mass loss will be the major reason for the deviation of the rising velocity at the lowest resolution in Fig. 5. Interface plot in the early stage of the rising process is shown in Fig. 7 with three different grid resolutions at simultaneous time. We can clearly see that the interface evolution is converging with increasing grid resolution. Considering the simulation results from Figs. 5-7, a grid resolution of 120 × 720 has been chosen as the working resolution for current numerical simulation tests. 3.2 Bubble rising simulation with top free surface The bubbles generated from the evaporation column will be collected in the steam drum and then separated from free surface. The separated steam should be exhausted without any droplets for the efficient operation of the HRSG system. The dominant forces determining the physical behavior of the interface motion would be the gravity, surface tension, and inertia. We have chosen initial diameter (D in Fig. 4), rising distance (H in Fig. 4), and ambient pressure as the major controlling parameters. The material properties of the bubble (vapor in Fig. 4) and water (liquid in Fig. 4) will be a function of the ambient pressure. Considering the general operational condition of the HRSG system, we have selected a few possible values of the ambient pressure for current parametric study. We have summarized the simulation condition for the various input pressure in Table 1. We have simulated a two-dimensional axisymmetric bubble rise problem with top free surface as described in Fig. 4.
5
10
20
ρL [kg/m ]
958.4
942.7
914.9
886.6
849.1
ρG [kg/m3]
0.5969
1.142
2.699
5.206
10.17
µL×104 [N-s/m]
2.819
2.308
1.795
1.498
1.257
µG×10 [N-s/m]
0.1227
0.1298
0.1407
0.1504
0.1617
σ [N/m]
0.05891
0.05484
0.04824
0.04208
0.03467
4
Fig. 7. Interface profile of the single bubble rise simulation for three difference grid resolutions.
2
Simulation condition is exactly the same as the single bubble rise case in section 3.1 except top free surface presence. The radial length of the simulation domain stays the same as single bubble rising simulation, but the height of the simulation domain changes according the rising distance, H. The properties of steam water at 5 atm, typical operating condition of a small HRSG system, has been chosen as the reference case. Property ratio becomes ρ* = 0.00295 for the density and µ* = 0.07839 for the viscosity, respectively. Initially, spherical bubble has been placed at (0, 3.5D) position. The box size of 3D × 11D has been used for the reference case. The resulting interface evolution has been plotted in Fig. 8. The rising bubble becomes elongated in longitudinal direction to form an ellipsoidal bubble, then hits the top free surface. The thin filamentary interface between bubble and top surface breaks up to form a satellite droplet at the top of the free surface. For the bubble rise with top surface, the focus has been to identify the satellite droplet formation after bubble merging to the top free surface. The formation and behavior of the satellite droplet after impact would be the critical factor for proper operation of the steam drum in the HRSG system. The best possible scenario would be no droplet generation since an entrained droplet at the upper steam collector would damage the system performance. The droplet might explode with oversaturation temperature or condense the surrounding steam to form larger droplet, which ultimately hinders steam flow. Even if there is a droplet formation after impact, sufficiently small size of the droplet could be entrained due to the magnitude of upward steam velocity from phase change process at the free surface. So the size of the satellite droplet would be the important output parameter determining the performance of steam drum system. In addition to the size of the satellite droplet, its velocity would also determine the destination of the droplet. With large velocity, the detached droplet would travel further to reach the steam collecting area. To analyze the behavior of the bubble with top free surface more relevantly for the HRSG system, we chose our output parameters as detaching time, time required for the rising bubble to interact with top free surface, diameter of detached droplet, maximum velocity of the detached droplet. We tested the grid convergence of the bubble rise simulation with top free surface. The non-dimensional grid size for the reference case (see Fig. 8) is 0.025, which has been identi-
3146
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
Fig. 8. Interface evolution of the single bubble rise with top surface at several time steps.
(a)
(b) Fig. 9. Grid convergence test for single bubble rise with top surface: (a) convergence of detaching time vs. grid size; (b) convergence of diameter of detached drop vs. grid size.
fied as the working resolution for the single bubble rise test. We further decreased or increased the grid size to investigate its effect on simulation results. The non-dimensional detaching time of the satellite droplet and non-dimensional diameter of the detached droplet vs. grid size is plotted in Fig. 9. As can be seen there, the working grid size of 0.025 is again appropriate for the bubble rise with top surface. The merging distance, dl, in Fig. 3 has been set to 0.025, single grid size, for
the all simulation considered. As discussed before, active controlling parameters for the droplet formation after merging of the rising bubble with top surface are rising distance, initial diameter of the bubble, and ambient pressure for the material properties. To identify the individual effect from each input variable of the rising distance, initial diameter of the bubble, and ambient pressure, we fixed the other two variables and varied only one variable of interest from the reference case. In Fig. 10, the effect of rising distance (H), ambient pressure (Patm), and initial bubble diameter (D) on non-dimensional diameter of detached droplet is shown. As can be seen, rising distance and ambient pressure have less effect than initial droplet on the size of detached droplet. With increasing distance, H, and ambient pressure, Patm, the size of the detached droplet is increasing a little but relatively small amount. The initial size of the bubble has the dominant effect on the size of the satellite droplet. Initial bubble, which is less than 1 mm or greater than 8 mm in diameter, doesn’t form any satellite droplet. It can be readily expected that there will be no detached droplet with small initial bubble. The inertial force from the small rising bubble is not enough to make a detached surface after impact. In Fig. 13, we have shown interface evolution with large initial diameter of the bubble (D = 30 mm). The rising bubble tends to elongate further with drag force in horizontal direction with increasing initial volume. So the lower portion of the bubble becomes hollow to make a spherical cap shaped interface. The bottom surface will ultimately hit the topside of the bubble to make a donut shaped interface as in Fig. 13. This ring type interface will rise further to merge with top free surface but does not have enough inertia to from a satellite droplet. It is interesting to note that there is a maximum of the detached droplet size in the middle of the considered size distribution of the initial bubble. The detaching time is exactly proportional to rising distance, H, as expected as described in Fig. 11. The bubble rising velocity is converging to a constant very quickly as can be seen from Fig. 5. The non-dimensional detaching time has very weak effect from the ambient pressure, similar to the size of the detached droplet. Again, initial diameter of the bubble is the major controlling parameter for the detaching time. The detaching time is increasing with larger initial bubble size. The maximum velocity of the detached droplet is depicted in Fig. 12 with each controlling parameter. Again, rising distance and ambient pressure have an insignificant effect on the maximum velocity of the detached droplet. With changing ambient pressure, there is a minimal value near the reference pressure of 5 atm even though relative changes are fairly low. Initial diameter of the bubble has a major effect on the final velocity of the detached droplet. In this case, there is a minimum value of the maximum velocity as opposed to the result of detached droplet size. It can be readily inferred from Fig. 10 that the size of the detached droplet is inversely proportional to the maximum velocity.
3147
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
(a)
(b)
(c)
Fig. 10. Effects of input variables on the non-dimensional diameter of detached drop: (a) effect of distance to the top surface; (b) effect of ambient pressure; (c) effect of initial bubble diameter.
(a)
(b)
(c)
Fig. 11. Effects of input variables on the non-dimensional detaching time: (a) effect of distance to the top surface; (b) effect of ambient pressure; (c) effect of initial bubble diameter.
(a)
(b)
(c)
Fig. 12. Effects of input variables on the non-dimensional maximum velocity of detached drop: (a) effect of distance to the top surface; (b) effect of ambient pressure; (c) effect of initial bubble diameter.
4. Conclusions The bubble rise problem with top free surface has been analyzed numerically with level contour reconstruction method. The LCRM is a hybrid type interface method for multiphase flows which combines some essential features of front tracking and level set techniques. By adopting high order tetramarching reconstruction, we can maintain excellent interface fidelity and perfectly interconnected interfacial elements with
complex topological changes. A benchmarking test of single bubble rising has been performed to show the accuracy of axisymmetric version of the LCRM. The results shows good convergence with the existing body fitted numerical solution. For the bubble rise with top free surface, the focus has been to identify the formation of the satellite droplet. The effect from the rising distance, ambient pressure, and initial diameter of the bubble on the size and velocity of the detached droplet and detaching time has been identified. The influence of the initial size of the bubble is dominant for the interface evolution with
3148
S. Shin / Journal of Mechanical Science and Technology 26 (10) (2012) 3141~3148
the top free surface. The size of the detached droplet has maximum value on specific initial bubble diameter as opposed to the maximum velocity of a detached droplet showing minimum value. The interface generally undergoes severe transition of the topological evolution. More quantitative research of multiple bubbles interaction with top free surface, especially in 3D, is necessary and currently underway.
Acknowledgment This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (20090071479, 20100016889).
References [1] G. Ryskin and L. G. Leal, Numerical simulation of freeboundary problems in fluid mechanics. Part 2, J. Fluid Mech., 148 (1984) 19-35. [2] A. Esmaeeli and G. Tryggvason, Direct numerical simulation of bubbly flows. Part 1. Low Reynolds number arrays, J. Fluid Mech., 377 (1998) 313-345. [3] R. Scardovelli and S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Ann. Rev. Fluid Mech., 31 (1999) 567-603. [4] S. Osher and R. P. Fedkiw, Level set methods: An overview and some recent results, J. Comput. Phys., 169 (2001) 463502. [5] H. G. Choi, A numerical study on the entry problem of an elastic body using a level set method and a monolithic formulation, Journal of Mechanical Science and Technology (2012) (in press). [6] H. G. Choi, A least square weighted residual method for level set formulation, Int. J. Numerical Methods in Fluids, 68 (2012) 887-904. [7] M. Sussman and E. G. Puckett, A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows, J. Comput. Phys., 162 (2000) 301-337. [8] G. Son and N. Hur, A coupled level set method and volume of fluid method for the buoyancy-driven motion of fluid particles, Int. J. of Numerical Heat Transfer, 42 (2002) 523-542. [9] J. Choi and G. Son, Numerical study of droplet motion in a microchannel with different contact angles, Journal of Mechanical Science and Technology, 22 (2008) 2590-2599. [10] D. Enright, R. Fedkiw, J. Ferziger and I. Mitchell, A hybrid particle level set method for improved interface capturing, J. Comput. Phys., 183 (2002) 83-116.
[11] E. Aulisa, S. Manservisi and R. Scardovelli, A mixed markers and volume-of-fluid method for the reconstruction and advection of interfaces in two-phase and free-boundary flows, J. Comput. Phys., 188 (2003) 611-639. [12] S. Shin and D. Juric, Modeling three-dimensional multiphase flow using a level contour reconstruction method for front tracking without connectivity, J. Comput. Phys., 180 (2002) 427-470. [13] S. Shin and D. Juric, High order level contour reconstruction method, Journal of Mechanical Science and Technology, 21 (2007) 311-326. [14] S. Shin and D. Juric, A hybrid interface method for threedimensional multiphase flows based on front tracking and level set technique, Int. J. Numerical Methods in Fluids, 60 (7) (2009) 753-778. [15] I. Yoon and S. Shin, Tetra-marching procedure for high order level contour reconstruction method, Advances in fluid mechanics VIII, WIT press (2010) 507-518. [16] J. Han and G. Tryggvason, Secondary breakup of axisymmetric liquid drops. I. Acceleration by a constant body force, Physics of Fluids, 11 (12) (1999) 3650-3667. [17] J. Chorin, Numerical solution of the Navier-Stokes equations, Math. Comput. 22 (1968) 745-762. [18] F. H. Harlow and J. E. Welch, Numerical calculation of time dependent viscous incompressible flow of fluid with free surface, Phys. Fluids. 8 (1965) 2182-2189. [19] C. W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes II, J. Comp. Phys. 83 (1989) 32-78. [20] M. Sussman, K. M. Smith, M. Y. Hussaini, M. Ohta and R. Zhi-Wei, A sharp interface method for incompressible twophase flows, J. Comp. Phys. 221 (2007) 469-505. [21] J. Hua, J. F. Stene and P. Lin, Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method, J. Comp. Phys. 227 (2008) 3358-3382.
Seungwon Shin received his B.S. and M.S. degrees in Mechanical Engineering from Seoul National University, Korea, in 1995 and 1998, respectively. He then received his Ph.D. from Georgia Tech in 2002. Dr. Shin is currently a Professor at the School of Mechanical and System Design Engineering at Hongik University in Seoul, Korea. His research interests include computational fluid dynamics, multiphase flow, surface tension effect, phase change process.
