MARCHING CUBES BASED FRONT TRACK ING METHOD AND ITS APPLICATION TO SOME INTERFACE INSTABILITY PROBLEMS*

WANG Jing-yi, ZOU Jian-feng, ZHENG Yao

Center for Engineering and Scientific Computation, School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China, E-mail: wangjingyi7219@126.com

REN An-lu

Department of Engineering Mechanics, School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China

MARCHING CUBES BASED FRONT TRACK ING METHOD AND ITS APPLICATION TO SOME INTERFACE INSTABILITY PROBLEMS*

WANG Jing-yi, ZOU Jian-feng, ZHENG Yao

Center for Engineering and Scientific Computation, School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China, E-mail: wangjingyi7219@126.com

REN An-lu

Department of Engineering Mechanics, School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China

A front trackingmethod based on amarching cubes isosurface extractor, which is related filter generating isosurfaces from a structured point set, is provided to achieve sharp resolution for the simulation of non-diffusive interfacial flow. Compared with the traditional topology processing procedure, the current front trackingmethod is easier to be implemented and presents high performance in terms of computational resources. The numerical tests for 2-D highly-shearing flows and 3-D bubblesmerging process are conducted to numerically examine the performance of the currentmethodology for tracking interfaces between two immiscible fluids. The Rayleigh-Taylor (RT) and Richtmyer-Meshkov (RM) instability problems are successfully investigated with the presentmarching cubes based front trackingmethod.

interfacial flow, front tracking,marching cubes, bubblesmerging process, Rayleigh-Taylor (RT) instability, Richtmyer-Meshkov (RM) instability

Introduction

The flow problems with two ormorematerials are frequently encountered in engineering and cause significant difficulties for numerical simulation. Examples of such flow problems include the Rayleigh-Taylor (RT) and Richtmyer-Meshkov (RM) instabilities generated and accelerated by body forces or shock waves. These instabilities are common in such applications as Inertial Confinement Fusion (ICF) and supernova explosions. The interfacial flow could be diffusive or non-diffusive, which is determined according to the spacial and temporal scales of the physical problem. For the diffusive case,mixing at themolecular level occurs in the vicinity of the interface, but the sharp interface of non-diffusive flow could remain in the whole course of front propagation. The characteristic temporal scale of the ICF and supernova explosions is 6-9 order shorter than themolecularmixing time. Hence, the present work concentrates on the numerical investigations of non-diffusive interfacial flow.

The current nonlinear theory has not been so powerful to analyze the physicalmechanism of the interfacial flow. Therefore Computational Fluid Dynamics (CFD) becomes themajor technique to understand the interfacial flow. The common numerical schemes cannot been applied directly to study nondiffusive interfacial flow: for simulation with a lower order scheme, serious numerical diffusion near the front w ill smear the sharp interface out, and for a higher order case, spurious perturbation w ill damage the interface geometry. For the past decades, a number of researchers have been attracted to study the special numerical techniques for the simulation of the interfacial flow. Such new techniques include the Volume Of Fluid (VOF), level set and front trackingmethods. The variants of VOFmethod developed in those years, such as Flux-Corrected Transport (FCT)[1], Least-Square Approach (LSA)[2], differing themselves in the reconstructing technique for interface geometry. Yang et al.[3]devised an Unsplit Lagrangian Advection (ULA) scheme for the VOFmethod, which is developed based on the algorithm of Piecew ise Linear Interface Construction (PLIC). The numerical results for three common tests (i.e., the shearing flow, RT instability and Zalesak rotation test) by using the ULA showed that the ULA scheme can achieve the accuracy higher than the SLA and FCT algorithms. By using the VOFmethod and its variants,many researchers have contributed to the investigation of complicated interfacial flow problems. The viscosity effects on the behavior of an initially spherical buoyancydriven bubble rising in an infinite and initially stationary liquid have been investigated by Wang et al.[4]. Yang et al.[5]developed amodel to study the highly nonlinear interactions between waves and rubblemound breakwaters. Zhang et al.[6]computed the flow and sedimentmovement in a cavity with the VOFmethod, capturing the time varying free surface. However the VOFmethod ismuchmore diffusive than the true physical process by intrinsic property termed from the VOF algorithm itself.

The level setmethod is another kind of volume trackingmethod, which defines the interface by a level set function (ψ). In thismethod, the topology processing is automatic, which is the excellent feature of the level setmethod. But several drawbacks limit its applications to a small range of interfacial flows. Itsmain drawbacks are: firstly, the level set function is not a physical quantity and nonphysical featuresmay occur in the interface interactions; and secondly, it is difficult to solve problems with three ormorematerials.

The front tracking approach is a numerical tool for the accurate resolution of sharp interface of the non-diffusive problems. Discontinued physical values of flow states are allowable in the two sides of front structures and no numerical diffusion w ill be introduced in the interface treating process.

Tryggvason et al.[7]have solved the interfacial flow under the context of incompressible fluid by the front trackingmethod. Themark function I( x), which is obtained from a special Poisson-type equation connecting the interface information and the field distribution, is used to identify fluids from each other. The introduction of I( x) functionmakes it possible to obtain a single field description of the flow with two ormorematerials and no interior boundaries are needed to handle in their solving stage, which is the same as the VOFmethod. Hence, this implementation of front trackingmethodmay have good convergence and stability features for flows with large density ratio. However on the other side of the problem, the usage of I( x) w ill create distinct numerical diffusion near the front and smooth the sharpness of the interface structures. Muradoglu and Tryggvason[8]continued to present a general computational framework based on the front-trackingmethod for simulations of interfacial flows with soluble surfactants. The front trackingmethod developed by Tryggvason et al.[8]has also been adopted to simulate the dynamic response of a 3-D bubbly cluster[9].

Based on compressible fluid equations, Glimm et al.[10]successfully implemented the front trackingmethod to solve the sharp contact discontinuity. Their implementations of front trackingmethod contained a large number of codes to handle themerging and detaching of triangular facets, which construct the complicated front structures. It is known that the topology processing takes themajor portion of the simulation time, and any enhancement on topology processing issue w ill be helpful to save the computational resources. A general-purpose front tracking package for the geometry and dynam ics of an interface has been developed and published by Du et al.[11]. They have also assessed the performance of this front tracking package, by comparing with publicly distributed interface codes (the level setmethod), and other published performance results (the VOF and othermethods).

Besides the research workmentioned above, Lallemand et al.[12], W itteveen et al.[13], Gois et al.[14]have also conducted distinguished improvements to the traditional front tracking approach and demonstrated its advantages in resolving discontinuities for interfacial flows. There has beenmany publications released to demonstrate the efficiency of fronttrackingmethod, for example, Terashima and Tryggvason[15]improved the front-trackingmethod, wheremarker points were used both for tracking fluid and constructing the Riemann problem on the interfaces, and coupled with ghost-fluidmethod to solve the interfacial state value. Muradoglu and Tasoglu[16]successfully simulated the impact and spreading of a viscous droplet on dry solid walls. Rawat and Zhong[17]resolved interactions of flow-disturbance with shock waves in a way with a new fixed grid setup for front tracking. Li et al.[18]used the front tracking code to simulate crystal grow th and solute precipitation, a Stefan problem which develops dynamically and requires high resolution of the interface geometry. Besides, the simulation ofmoving interfaces such as a fire front poses a challenge due to the fact that a large scales of space and time often require small time steps, Filippi et al.[19]presented a novelmethod based on the front trackingmethod, thus enables the simulation of large-scale fire-spread systems. There aremany other new publications further demonstrated the capability of the efficiency of fronttrackingmethod.

The volume trackingmethods solve a set of single field equations and is easy to be implemented,but this type ofmethods are diffusive and the sharpness of the captured interfaces is poor. The front trackingmethods are non-diffusive, which are in accord with the real physical process, and are the prom ising tools to achieve accurate and sharp interface geometry.

Here, we choose the front trackingmethod as our tool to precisely investigate the complex surface instability problem. The original topology processing procedure has beenmodified in our work tomeet the requirement of high performance simulation. In our simulation, the fluids are indexed with two ormoremarker components (0, 1...), and themarching cubes algorithm is adopted to compute the interface from the discrete field ofmarker components. Hence, the topology processing in our simulation is conducted implicitly and no additional operations are needed to consider themerging and detaching of fluid fronts. A detailed description of the front propagation procedure for one whole time step is given in Section 2. Section 1 presents themathematical form of the front trackingmethod. Section 2 presents themarching cubes based topology processing procedure in detail. The simulation results on the RT and RMsurface instability problems are given in Section 3. And finally, Section 4 is a brief conclusion of our research work.

1. Mathematical form of front trackingmethod

The Navier-Stokes equations can be w ritten in a conservation form as follows

Here τ is the viscous stress tensor, which is related to the rate of strain by the constitutive equation in the form of components

And the viscous coefficient μ can be calculated using Sutherland’s law, λ=?2/3μ, and the heat flux terms are given by Fourier’s law as q= (qx, qy, qz)=?k? T, with k being the heat transfer coefficient and T the flow temperature.

The front trackingmethod is a promising approach to achieve sharp resolution of interfaces between distinctmaterials. In our implementation of the front trackingmethod two sets ofmeshes are defined: the 3-D volumemeshes for solving flow equations and the 2-D surfacemeshes (i.e., the triangular facets constructing the interface geometry) formodeling the interface state propagation. The fluid interfaces divide the full flow domain into two ormore sub-domains and the interface segments act as the interior boundaries when the flow equations are solved in each sub-domain. The propagation of the 2-D surfacemesh points obey the Rankine-Hugoniot rule. In the course of front propagation, the interfacetopology w ill change (merge or detach) and a certain operation should be adopted to handle it to avoid nonphysical interface segments. The topology processing procedure plays an important role in the implementation of front trackingmethod. The presentmarchingcubes-based topology processing procedure is given in Section 2 in detail.

2. Marching-cubes based topology processing procedure

2.1 Main steps to treat interfaces

After a certain time steps, the line segments or triangular facets on the fluid front propagate according to the local fluid states and tend to tangle with each other. Hence the topology processing operations should be performed frequently to avoid nonphysical fluid states. In our front tracking implementation, the complete interface treating process includes threemain steps: intersection computation,marker component detection and isosurface extraction.

(1) Intersection computation

Firstly the positions at which the triangular facets and the volumemeshes intersect should be computed. Figure 1 shows the 2-D case at the end of several steps of propagation, where one rising bubble or drop penetrates the other. The small color triangles in Fig.1 represent the intersections and the colors (black and white in Fig.1) tell us the exactmaterial components on both sides of the front.

Fig.1 Intersection computation for 2-D case

Fig.2 Marker component detection for 2-D case

(2) Marker components detection

Here, themarker component for eachmesh point is determ ined by the previously computed intersections. In our code, the depth-first-searching algorithm is used to loop over all themesh edges and achieve the right fluid component field. For themesh element containing no interface segments, the vertices of each edge share the samemarker component with external boundary (see themesh edge identified as IAJA→IBJAin Fig.2). For the intersecting edge with untangled interface segments (IBJC→ICJCin Fig.2), the components of the two vertex points of this edge can be determ ined from the color triangles. For the case of tangled interface (IEJC→IEJDin Fig.2),more care should be taken to detect and delete nonphysical intersecting points. When a full searching loop is finished, the discrete component field has become consistent. In Fig.2, the fluid component of eachmesh point is indicated by a color circle.

(3) Isosurface extraction

Till now, the full flow field has consistent components and it is ready to construct new interface structure from the distribution of fluid components.

The fluid front is thought of the isosurface of the discretemarker components (0,1...) and a robust isosurface extractor contained in the Visualization Toolkit (VTK) package[20]is applied here. The VTK package is a C++ open source package for the computer visualization, whichmakes the isosurface generating procedure easy to be coded, and the complicated topological treating process is avoided in our front tracking operation. As our 3-D volumemeshes are composed of structured points, themarching cubes implementation in the VTK package, that is, the C++ class VTKmarching cubes, is selected as the isosurface extractor, which is themost frequently referenced isosurface extractor. To achieve the isosurface at the new time step, we shouldmake a loop over all rectangular (for 2-D case) or hexahedral (for 3-D case) elements, and generate the local topology in each one. Care should been taken to resolve the ambiguous configurations, when the local topology segments are composed to generate the global fluid front, which occurs when the front segment in one element cannot be uniquely determined and the neighboring elements should be considered in a coupledmanner. The asymptotic decider proposed by Nielson[21]has solved the ambiguous configurations very well and w ill be adopted in our work.

By the use of VTKmarching cubes object in the VTK package, the isosurface with value 0.5 can be achieved (see Fig.3(a) for 2-D case). The achieved isosurface of value 0.5 and themesh element intersect at them idpoint of themesh edge. Therefore, the previous computed intersecting information in the “Intersection computation” step should be used tomodifythe isosurface produced by the VTK. The final updated front structure is given in Fig.3(b).

Fig.3 Isosurface extraction for 2-D case

By using themarching-cubes-based interface treating procedurementioned above, the fluid front is updated and the interfacial flow simulationmarches to another time step.

2.2 Numerical tests for 2-D and 3-Dmoving interfaces

In this sub-section 2-D and 3-D tests are conducted to assess the currentmethodology for tracking interfaces between two imm iscible fluids. The 2-D test is designed to investigate the accuracy of a sharp interface involved in a complicated shear flow. And the accuracy of interfacemerging procedures is validated by a 3-D bubble-merging computation.

Fig.4 Interface evolution in a highly-shearing flow field

The first numerical test for the tracking procedure is performed in a highly-shearing flow field. In this case, the velocity is prescribed as u=αsin(πx)· sin(2πy) and v=βsin(π y)sin(π y)sin(2π x), where x∈[0m,1m], y∈[0m,1m]. The interface is firstly composed of four bubblesmoving at the prescribed velocity (α=?1and β=1) for one second. And then the velocities are reversed for another second (α=1 and β=?1).

The series of results of this 2-D test is presented in Fig.4 with the time increment of 0.1 s. The bubbles aremoving under the influence of fluid shear and perform considerable deformation. The current front tracking approach has resulted in the same initial interface at the time t=2.0 s, despite being highly distorted at t=1.0 s.

Fig.5 3-D interfacemerging process

Another 3-D test is designed to investigate themergingmotion of four bubbles whichmove according to a velocity defined by their own normal vector. Figure 5 presents the time history of themerging process in detail.

3. Simulation results of flow instability problems

3.1 Single-mode RT instability simulation

The RT instability occurs when a light fluid is pushed by a heavy one. The perturbations evolve in a strong nonlinearmanner, which play an important role in the investigation of ICF problems. Here, we apply themarching-cubes-based front trackingmethod to simulate the RT instability with single andmultiplemodes.

A single-mode RT simulation was conducted to validate the code ofmarching-cubes-based interface treating procedure. Consider a 2-D problem with the computational lim its in the x and y directions are 0.01m and 0.1m, respectively. The left and right boundaries are set to be periodic. The initial disturbance has the follow ing form

where W is the box w idth.

Amesh with 100 points in the x direction and 1 000 in the y direction is used in the computation.The density ratio of the light to the heavy fluid is 1 to 2 (the light fluid is above the heavy one) and the corresponding Atwood number is 0.333, which is defined as

wherehρ denotes density of the heavy fluid andlρ that of the light one. The gravitational acceleration is set to be 0.14m/s2.

Fig.6 Velocity distributions at t =124.5ms and t= 249ms

Fig.7 Time history of fluid front structures

Figure 6(a) presents the symmetrical velocity distributions in the early phase of time 124.5ms. The rising bubbles and falling spikes are well-shaped. But as the timemarches to 249ms, the flows are fully unsteady and become turbulent (see Fig.6(b)). Figure 7 gives the fluid fronts on the time slices of 167ms, 188ms, 210ms, 230ms and 249ms, respectively. The unsteady roll-ups become distinct at the time of 210ms as shown in the third plot of Fig.7, not only in the rising bubbles, but also in the falling spikes. The formation, rollup and pairings of the interfaces have been clearly captured along the fluid front by using the present front trackingmethod, and even the smallscale vortices are well resolved.

The expression h=αAgt2is deduced from the linear theory and used to define the initial RTmixing for the height h of the rising bubble. Themixing rate α has the experimental values of 0.063± 0.013 (refer to [22]). The amplitudes of bubbles and spikes are plotted as the function of the time in Fig.8. In our simulation, the value of α is taken as 0.078, which is very close to the experimental data.

Fig.8 Amplitudes of bubbles and drops vs time

3.2 Multi-mode RT instability simulation

In this simulation, the computational domain is a 0.02m×0.02m×0.08m volume, and there are 100, 100 and 400 intervals in its three directions, respectively. The density ratio of the light to the heavy fluid is 1 to 3 (the light fluid is above the heavy one). Themulti-mode simulations are initialized with random surface perturbation, which could be expressed with the formulation of Eq.(7), to investigate the RT instability dominated bymode coupling. Themodes containing resolvable wave numbers between 16 π and 32 π are assigned with random finite amplitudes. A ll simulations are parallelized with the MPI library on our SGImachine.

The 3-D interface structures at time 55ms, 95ms, 135ms and 170ms are given in Fig.9. The strongmode coupling features can be found in the later phase of the interface evolution.

Figure 10 shows the surfacemeshes of the fluid front at time 135ms. The presentmarching-cubesbased front trackingmethod solves the topological change automatically, and sufficientmesh points are guaranteed on the fluid front to capture themulti-scale interface structures.

Fig.9 Fluid front evolution for the 3-Dmulti-mode RT problem

Fig.10 Lower surfacemesh on the front structure

3.3 RMinstability simulation

The purpose of the third simulation is to investigate the RMinstability induced by an incident shock wavemoving across the interface between two different fluids. In our case, the computational domain is a 0.01m×0.04m rectangular area, and there are 100 and 400 intervals in its two directions, respectively. The density ratio of the light to the heavy fluid is 1:2 (the light fluid is below the heavy one, and its density is 1 kg/m3). Suppose that the upper side of the domain is the wall of the shock tube, and the other side is directly linked to open air. The calculations are performed in Cartesian coordinates. At t=0μ s, the shock wave is at the position of y=0.01m , and ismoving upward with the Mach number of 2.0. The static pressure in both fluids is 1 atm, and γ in both fluids is 1.67. The initial interface between the two fluids is given by

As Fig.11(a) shows, the white component is lighter one while the black component is heavier one. The top wall is the Neumann boundary, the bottom wall is the Dirichlet boundary and the two sides of the domain are periodic. Figure 11(b) shows themoment before the incident shock wavemoves across the interface, and then the shock wave wouldmake interface amplitude become larger. Figure 11(c) shows themoment before the shock wave is about tomove across the interface after the shock wave re-bounced back from the top wall. Figures 11(d) and 11(e), after being re-shocked, the second shock produces a light gas incident on a heavy gas configuration and an inversion of the interface occurs.

To perform a further detailed research on the RMinstability, the front trackingmethod was applied to simulate a shock-tube experiment which was conducted by Jacob and Krivets[23]. In Jacob and Krivets’s experiment, the shock tube is a square test section with a cross section of 0.089m×0.089m×0.075m. Air and SF6were used as work fluids. In our simulation, 2-D assumption was applied and the computational domain was simplified as a rectangular of 0.06m×0.9m. We keep the initial amplitude and wavelength the same as in the experiment, as well as other parameters.

The present perturbation is sinusoidal with wave number k and its amplitude grow th rate in the RMinstability can be approximately expressed as ν0= kΔ uA+, where Δu is the change in translational velocity of the interface due to the action of the shock, A+=(?)/(+)is the post-shock Atwood number, and a0+is post-shock initial perturbation amplitude. In the present simulated case, Δu= 68.6m/s. The computational parameters of two work fluids, air and SF6, are presented in Table 1.

The simulation had a grid resolution of 260× 3 900 corresponding to the actual experimental domain of 0.06m×0.9m and its CFL number is 0.75. Because the length of the computed rectangular is longer than that in the experiment, we are only going to study the interface evolving before the re-shock onto the interface, so there would be no phase inversion where the re-shock wave is reflected backward from the upper wall.

In Fig.12, the time histories start with themoment the incident shock impinge onto the interface. The fluidmarked with black component is SF6andthe white component is air. The small-scale flow structures like roll-ups and pairings are caused by themisalignment of the pressure and density gradients across the contact between two differentmaterials[24].

Fig.11 Time history of front structures and shock wave

Table 1 Com putational parameters

Fig.12 Time histories of simulated component field at selected times

The simulation amplitude is compared to the experimental amplitude data in Fig.13. Though there is slight difference between the simulation and experimental data, the simulation data still keep a relatively consistent trend with the experimental data before reshock. The subsequent discrepancy is due to the arrival of the driver-based expansion wave notmodeled in the simulation.

Fig.13 Comparison of simulation and experimental amplitude data

4. Conclusion

The front trackingmethod based onmarching cubes has been developed and validated in the present work. We have first implemented the front tracking approach to simulate one 2-D numerical test (interface evolution in a highly-shearing flow field) and one 3-D (multiple bubblesmerging) numerical test to demonstrate the accuracy and effectiveness of the front trackingmethod for interface tracking problem. Then we have a brief numerical study for the RT instability and obtain a satisfactory result for the evolution of RT instability. The simulation of themulti-mode instability problem, initialized by random perturbation, exhibits a detailed evolutionmechanism of coupling bubbles and drops. Ourmain goal is to use this present front tracking code to simulate the development of RMinstability, the small-scale flow structures like roll-ups and pairings are well simulated and the numerical and experimental data are in a good agreement up to the time relatively late but before re-shock, though the amplitude between the spike and bubble from simulation is slightly larger than the amplitude inexperiment. Note that in the experiment a weak expansion wave is generated by the incident shock when it passes over the slots, and this interaction subsequently results in a reduction in the interface velocity and strength of the transm itted shock wave, whichmainly lead to the difference between simulation and experiment. In order to obtain amore detailed analysis for the RMinstability, we would further utilize the present front tracking code to simulate the development of this phenomenon and combine the experiment and simulation data with existing predictionmodel for the RMinstability.

Acknowledgements

We also appreciate the helpful discussions among themembers of the CFD group at the Center for Engineering and Scientific Computation, Zhejiang University.

[1] RUNDMANMM. Volume-trackingmethods for interfacial flow calculations[J]. International Journal for Numerical Methods in Fluids, 1997, 24(7): 671-691.

[2] WACHEMB. G. M., SCHOUTEN J. C. Experimental validation of 3-D Lagrangian VOFmodel: Bubble shape and rise velocity[J]. AIChE, 2002, 48(12): 2744- 2753.

[3] YANG Wei, LIU Shu-hong and WU Yu-lin. An unsplit lagrangian advection scheme for volume of fluidmethod[J]. Journal of Hydrodynam ics, 2010, 22(1): 73-80.

[4] WANG Han, ZHANG Zhen-yu and YANG Yong-m ing et al. Viscosity effects on the behavior of a rising bubble[J]. Journal of Hydrodynam ics, 2010, 22(1): 81-89.

[5] YANG Chi, LU Hai-dong and L?HNER Rainald. On the simulation of highly nonlinear wave-breakwater interactions[J]. Journal of Hydrohynam ics, 2010, 22(5 Suppl.): 975-981.

[6] ZHANG Xiao-feng, CUI Zhan-feng and LU Xin-hua. Simulation of sediment deposition in a cavity with free surface[J]. Journal of Hydrodynam ics, 2010, 22(5): 626-633.

[7] TRYGGVASON G., BUNNER B. and ESMAEELI A. et al. A front-trackingmethod for the computations ofmultiphase flow[J]. Journal of Com putational Physics, 2001, 169(2): 708-759.

[8] MURADOGLU M., TRYGGVASON G. A. front-trackingmethod for computation of interfacial flows with soluble surfactants[J]. Journal of Computational Physics, 2008, 227(4): 2238-2262.

[9] CHEN Hui, LI Sheng-cai and ZUO Zhi-gang et al. Direct numerical simulation of bubble-cluster’s dynamic characteristics[J]. Journal of Hydrodynam ics, 2008, 20(6): 689-695.

[10] GLIMMJ., GROVE J. W. and LI X. L. et al. Threedimensional front tracking[J]. SIAMJournal of Sci- entific Computing, 1998, 19(3): 703-727.

[11] DU J., FIX B. and GLIMMJ. et al. A simple package for front tracking[J]. Journal of Com putational Physics, 2006, 213(2): 613-628.

[12] LALLEMAND P., LUO L. S. and PENG Y. A lattice Boltzmann front-trackingmethod for interface dynam ics with surface tension in two dimensions[J]. Journal of Com putational Physics, 2007, 226(2): 1367-1384.

[13] WITTEVEEN J. A. S., KOREN B. and BAKKER P. G. An improved front trackingmethod for the Euler equations[J]. Journal of Computational Physics, 2007, 224(2): 712-728.

[14] GOIS J. P., NAKANO A. and NONATO L. G. et al. Front tracking withmoving-least-squares surfaces[J]. Journal of Com putational Physics, 2008, 227(22): 9643-9669.

[15] TERASHIMA H., TRYGGVASON G. A front-trackingmethod with projected interface conditions for compressiblemulti-fluid flows[J]. Computers and fluids, 2010, 39(10): 1804-1814.

[16] MURADOGLU M., TASOGLU S. A front-trackingmethod for computationalmodeling of impact and spreading of viscous droplets on solid walls[J]. Computers and Fluids, 2010, 39(4): 615-625.

[17] RAWAT P. S., ZHONG X. L. On high-order shockfitting and front-tracking schemes for numerical simulation of shock-disturbance interactions[J]. Journal of Com putational Physics, 2010, 229(19): 6744-6780.

[18] LI X. L., GLIMMJ. and JIAO X. M. et al. Study of crystal grow th and solute precipitation trough front trackingmethod[J]. Acta Mathematica Scientia, 2010, 30B(2): 377-390.

[19] FILIPPI J. B., MORANDINI F. and BALBI J. H. et al. Discrete event front-tracking simulation of a physical fire-spreadmodel[J]. Journal of Simulation, 2010, 86(10): 629-646.

[20] SCHROEDER W., MARTIN K. and LORENSEN B. The visualization toolkit: An object-oriented approach To 3D graphics[M]. New York: Kitware, Inc., 2003.

[21] NIELSON G. M., HAMANN B. The asymptotic decider: Resolving the ambiguity inmarching cubes[C]. Proceedings of Visualization’91. San Diego, CA, USA, 1991, 83-91.

[22] YOUNGS D. L. Turbulentmixing due to Rayleigh-Taylor instability[C]. American Physical Society Division of Fluid Dynam ics Meeting. San Antonio, TX, USA, 2008.

[23] JACOBS J. W., KRIVETS V. V. Experiments on the late-time development of single-mode Richtmyer-Meshkov instability[J]. Physics of Fluids, 2005, 17(3): 034105.

[24] ZABUSKY N. J. Vortex paradigm for accelerated inhomogeneous flows: Visiometrics for the Rayleigh-Taylor and Richtmyer-Meshkov environments[J]. Annual Review of Fluid Mechanics, 1999, 31(1): 495-536.

January 23, 2011, Revised May 13, 2011)

10.1016/S1001-6058(10)60152-7

* Project supported by the National Natural Science Foundation of China (Grant No. 10702064).

Biography: WANG Jing-yi (1987-), Female, Master Candidate

ZOU Jian-Feng, E-mail: zoujianfeng@zju.edu.cn