Interface reconstruction and advection schemes for volume of fluid method in axisymmetric coordinates
II NTERFACE R ECONSTRUCTION AND A DVECTION SCHEMES FOR V OLUME OF FLUID METHOD IN A XISYMMETRIC COORDINATES
A P
REPRINT
Ananthan M
Department of Mechanical EngineeringIndian Institute of ScienceBangalore, India [email protected]
Gaurav Tomar ∗ Department of Mechanical EngineeringIndian Institute of ScienceBangalore, India [email protected]
February 2, 2021 A BSTRACT
Volume of fluid(VOF) method is a sharp interface method employed for simulations of two phase flows.Interface in VOF is usually represented using piecewise linear line segments in each computationalgrid based on the volume fraction field. While VOF for cartesian coordinates conserve mass exactly,existing algorithms do not show machine-precision mass conservation for axisymmetric coordinatesystems. In this work, we propose analytic formulae for interface reconstruction in axisymmetriccoordinates, similar to those proposed by Scardovelli and Zaleski (J. Comput. Phys. 2000) forcartesian coordinates. We also propose modifications to the existing advection schemes in VOF foraxisymmetric coordinates to obtain higher accuracy in mass conservation. K eywords volume of fluid · axisymmetric coordinates · multiphase simulation Multiphase flows are ubiquitious in several industrial applications. In the last three decades, there has been a surge inthe numerical methods and algorithms for simulations of complex multiphase flows. There have been several differenttypes of interface capturing strategies that have been proposed for two-phase flows. The most popular of these are theLevel set method, Volume of Fluid (VOF) method, and Front tracking scheme [1, 2]. VOF methods with geometricadvection strictly conserve the volume of the two phases.Several improvements have been made since the inception of the method (see Hirt and Nichols[3, 4, 5, 6, 7, 8]). Themost simplest and earliest representation for interface reconstruction is simple line interface calculation(SLIC) in whichthe interface is approximated by horizontal or vertical lines. Subsequently, piecewise line interface construction (PLIC)was introduced where the interface is approximated as a linear line at an angle in the cell [9]. Higher order interfaceconstruction have been proposed (such as Parabolic reconstruction by [4]), but considering the associated computationalcost and complexity for geometric advection, PLIC is usually preferred.Scradovilli and Zaleski[10] proposed analytical formulae for the piecewise linear reconstruction of the interface incartesian coordinates that led to a significant speedup over the earlier iterative schemes. However, for curvilinearcoordinates (such as axisymmetric coordinate system), the proposed analytical formulae cannot be employed directly.In the present work, we derive similar analytic formulae for axisymmetric coordinates, which result in a speedup of ∼ over the iterative counterparts (Brent’s root finding method). Further, we demonstrate that the existing interfaceadvection schemes in VOF for axisymmetric coordinates are not strictly mass conserving. In this study, we proposemodifications to the current operator split algorithms that result in machine-precision mass conservation in axisymmetriccoordinates. We show the efficacy of the proposed algorithms using several test cases. ∗ Corresponding author. a r X i v : . [ phy s i c s . c o m p - ph ] F e b PREPRINT - F
EBRUARY
2, 2021The paper is organized as follows. We first present analytical formulae for the interface reconstruction schemes inaxisymmetric coordinates in section 2. In section 3, we propose modifications in the existing interface advectionalgorithm for axisymmetric VOF and present test cases to show the efficacy of the scheme. Finally, in section 4, wediscuss the important conclusions.
Interface reconstruction in the volume of fluid(VOF) method requires the volume fraction field. Using the volumefraction field, a piece-wise linear or a higher order interface is constructed in a given grid-cell. Interface reconstructionis an integral part of the geometric advection schemes to ensure mass conservation property of the VOF method [2].Initial condition for a multiphase flow simulation requires the initial distribution of the volume fraction field, usuallyprovided as an implicit function of the spatial coordinates. VOFI library[11] is an open source library to initialise theliquid volume fraction field in cartesian coordinate systems accurately. In VOFI, for cells cut by the interface (see figure1), PLIC reconstruction method [9, 6] is employed to approximate the interface as a linear line segment, m . x = a, (1)where m is the local normal at the interface, x is a point on the plane, and a is the normal distance of the origin fromthe plane. Analytical relation, given by Scardovelli and Zaleski[10], between the volume fraction and the line constantis employed to determine the line constant a . Thus, for two-dimensional and three-dimensional Cartesian coordinatesystems, VOFI library can be directly employed for accurate assignment of the initial volume fraction field on a givendiscretized domain using an implicit equation of the interface. However, in curvilinear coordinate systems, for a givenimplicit function, the piece-wise linear interface constructed from VOFI would require computation of the volumefraction field using a formula specific to the curvilinear coordinates. For instance, for axisymmetric coordinate system,the modified Gauss area (shoelace) formula for computation of the area of a convex polygon is given by, V = π (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n (cid:88) i =1 ( x i + x i +1 ) ( x i y i +1 − x i +1 y i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (2)where ( x i , y i ) for i = 1 , ..., n (with x n +1 = x and y n +1 = y ) are the coordinates of the vertices of a convex polygonordered counter clockwise as shown in the figure 1.Figure 1: Coordinates of the vertices of a simple polygon cut by the interface ordered counter clockwise. Here x , y and x , y represents the end points of the PLIC line segment with interface normal m and line constant a .Thus, for initialization of the volume fraction field, C , once we obtain the linear interface in each grid cell using theVOFI library, we use the above formula to compute the volume, V , and assign volume fraction in each grid cell as, C = V πr c ∆ x ∆ y . (3)Here r c is the distance of the center of the cell from the axis of symmetry, and ∆ x and ∆ y are the grid-cell sizes in theradial ( r ) and axial ( y ) directions, respectively. We note that the above procedure is followed essentially to minimizethe error in the volume-fraction during initialization. To illustrate this, we initialise a torus of minor radius r t = 0 . PREPRINT - F
EBRUARY
2, 2021Table 1: Results for relative error in volume during initialization of a torus of radius, r t = 0 . , for different grid sizes.Relative error in volume: E = | V − V t | V t Grid Current Solver Gerris Solver ×
16 5 . × − . × − ×
32 7 . × − . × − ×
64 1 . × − . × − ×
128 5 . × − . × − and major radius r = 0 . in the center of a computational domain of size × as shown in the figure 7. The volumeof the torus can be analytically computed as V t = 2 π rr t , where the major radius, r , is the distance to the center ofthe torus from the axis of symmetry. We compare the results for various grid sizes with the results obtained using thepopular VOF based open source flow solver, Gerris[12], given in table 1. r y Figure 2: Reconstructed interface for 16 mesh points of a torus of radius r t = 0 . and major axis r = 0 . initialisedat the center of × domain where the dots represent the mid–points of the reconstructed PLIC line segments.Thus, we have shown that for curvilinear coordinates volume fraction field can be initialized up to machine accuracy.Now, we derive an analytic relation for PLIC reconstruction for axisymmetric coordinates on the lines of Scardovelliand Zaleski [10]. We use Youngs method[9] to get the interface normal ( m in equation.1) from fluid-1 ( C = 1 ) tofluid-2 C = 0 and is given by m = −∇ C/ |∇ C | . To complete the PLIC interface reconstruction for a given C , inaddition to the normal m , we also need to obtain the line constant, a , which is the normal distance of the interface fromone of the vertices of the computational cell. In what follows, we present a methodology to get the line constant( a )analytically for a given interface normal vector and the volume fraction of a mixed cell.As discussed in [10], using an analytical relation between the volume fraction( C ), interface normal ( m ) and the lineconstant ( a ), we can implement an if − else − if − end − if construct to determine the line constant, a . This approachis computationally much more efficient compared to the alternative iterative approach to get the line constant. Giventhe equation of the interface, m x + m y = a , all combinations of m , m (such that m + m = 1 ) can be reducedto one of the cases shown in figure 3, either by changing the origin or by changing the reference fluid from fluid- to fluid- , such that both m and m are positive and the left bottom corner of the mixed cell under consideration iscontained in fluid- . Figure 3 shows all the possible configurations for interface arrangement with m ≥ and m ≥ .We first discuss Case A shown in the figure 3 and similar procedure can be followed to obtain relations for other cases.For the axisymmetric coordinate system, the shaded area shown in the figure 3 for Case A is given by, V = π x + x ) m m − π x + x ) m m − πx y . (4)Using the equation of line m x + m y = a we have x = a/m − ( m ∆ y ) /m and x = a/m . In the presentstudy, we assume ∆ x = ∆ y , but the same analysis can be easily extended for ∆ x (cid:54) = ∆ y . Substituting x and x in the3 PREPRINT - F
EBRUARY
2, 2021Figure 3: Various cases which can arise in the standard configuration of the interface with m , m ≥ and fluid isoccupying the bottom left corner of the cell.equation 4 and collecting the terms in powers of a , we obtain, (cid:18) π ∆ ym (cid:19) a + (cid:18) − πm ∆ y m + 2 π ∆ yx m (cid:19) a + (cid:18) πm ∆ y m − π ∆ y m x m (cid:19) = V. (5)Thus, we have an analytical relation between the volume and the line constant a . We note that the above relationholds true only when the interface cuts through the top and the bottom edges of the cell shown for Case A in Figure 3: ( x − x ) ≤ ∆ x , y = ∆ y and y = 0 . These conditions yield the bounds on the values of a : m ∆ y ≤ a ≤ m ∆ x .Substituting the above bounds for a in the equation 5, yield the bounds on the limiting volumes: V = π ∆ y (cid:0) x m + 3∆ xm (2 m x − ∆ ym ) + ∆ ym (∆ ym − m x ) (cid:1) m . (6)and V = π ∆ y m (∆ ym + 3 m x )3 m (7)For a given volume fraction C and interface normal ( m , m ), the volume occupied by fluid- in the configurationsshown in the Fig.3 is given by: V = 2 πr ∆ x ∆ yC where r = x + ∆ x/ is the distance from the axis of symmetryto the cell center. If V ≤ V ≤ V then the analytical relation given by equation 5 can be used to determine the lineconstant a . For the quadratic equation in a given by equation. 5, we note that only one of the roots will satisfy therequired bounds on a for case A. For cases B and C, we obtain cubic equations that can be solved for a using Cardano’sformula or using Brent’s method to find the appropriate root with necessary bounds for the line constant. Following thesame approach we can get bounds on volume for case D as: V = π ∆ y ( − ym + 3∆ ym − m x + 6 m x )3 m (8)and V = π ∆ y m (∆ y + 3 x )3 m . (9)Figure 4 shows all the possible configurations of the interface and the volume bounds which separate each case as theinterface normal in radial direction varies from minimum to a maximum. Figure 4 clearly shows that the various bounds4 PREPRINT - F
EBRUARY
2, 2021 m V Case D Case ACase BCase CV V V V Figure 4: Bounds for different cases shown in figure 3 as a function of the x-component of the interface normal, m .The interfacial cell is placed at a distance of from the axis of symmetry with ∆ x = ∆ y = 1 . The total volume of thecell is V cell = 2 π given by the top boundary in the plot.for the cases shown in figure 3 do not overlap and provide a unique criterion for computing the line constant a . Thus,we can use the following algorithm to classify each case. Data:
V, m Result:
Identification of the case in which the standard interface belongs to. if m ≥ √ thenif V ≥ V then Case B else if V ≤ V then Case C else
Case A // V > V > V endelseif V ≥ V then Case B else if V ≤ V then Case C else
Case D // V > V > V endend Algorithm 1: Classification of the standard case of the reconstructed interface where m , m ≥ .We list below the analytical relation between the line constant( a ) and the volume( V ) for each case: Case A (cid:18) π ∆ ym (cid:19) a + (cid:18) − πm ∆ y m + 2 π ∆ yx m (cid:19) a + (cid:18) πm ∆ y m − π ∆ y m x m (cid:19) = V. (10)5 PREPRINT - F
EBRUARY
2, 2021
Case B − (cid:18) π m m (cid:19) a + (cid:18) π ∆ ym − πx m m (cid:19) a + (cid:32) π (∆ y + x ) m − πm ∆ y m + 2 π ∆ yx m − πx m (cid:33) a + (cid:32) − π ∆ y (∆ y + x ) m m + πm ∆ y m − π ∆ y m x m − πm x m + πm (∆ y + x ) m (cid:33) = V (11) Case C (cid:18) π m m (cid:19) a + (cid:18) πx m m (cid:19) a = V (12) Case
D a = 2 π ∆ y m + 3 π ∆ y m x + 3 m V π ∆ y (∆ y + 2 x ) (13)We note here that the other cases can be readily transformed into one of the cases listed in the figure 3 by either changingthe fluid (by using ( − C ) instead of C to compute the volume and inverting the interface normal m ) or by changingthe origin (keeping the location of the axis-of-symmetry same but inverting its direction).We now compare the above described analytical method with the iterative method for finding the line constant for thecase A given in figure 3. The relative error in the line constant is given in table 2 for the analytical and iterative methodswith different tolerances. We note that the iterative method to reach an error with a tolerance of − is about slowercompared to the analytical method.Table 2: Relative error in line constant a and the ratio of CPU time required by iterative method(Brent’s algorithm) tothat required by analytical method for repetitions for case A shown in figure 3 for different tolerances used in theiterative method. Comparision between the analytical and iterative reconstruction methodsMethod Tolerance Relative Error t iterative /t analytical Analytical − Iterative . × − . × − . Iterative . × − . × − . Iterative . × − . × − . Once the line constant, a , is obtained, the position of the endpoints of the linear approximation of the interface can becomputed thus completing the construction of a planar interface in a given computational cell. As discussed earlier, thismore precise description of the interface within the grid cell allows geometric advection which gives the VOF methodits strict mass conservation property while maintaining a sharp interface.In what follows, we discuss an operator split algorithm for the geometric advection of the interface in axisymmetriccoordinates. We note that the straightforward extension of the 2D cartesian algorithm does not yield accurate results, asalso indicated by the results obtained from the existing open source codes. We present here a scheme for accurate geometric advection of the volume fraction in the axisymmetric coordinates.We have used a uniform grid to describe the variables with volume fraction being stored at the cell centers ( C i,j ). Theincompressible fluid flow is determined by the velocity field which is defined at the cell faces ( u i +1 / ,j , v i,j +1 / ). Here, u denotes the radial direction velocity and v is the axial velocity. The velocity field satisfies the discrete divergence freecondition given by: ( ru ) i + ,j − ( ru ) i − ,j r i ∆ x + v i,j + − v i,j − ∆ y = 0 . (14)6 PREPRINT - F
EBRUARY
2, 2021Motion of the interface is governed by the advection equation for the volume fraction field, ∂C∂t + u · ∇ C = 0 (15)For incompressible fluids, conservation of the individual volumes of the two fluids results in the conservation ofmass. Thus, in the volume of fluid method, geometric advection of the volume fraction field is expected to yieldmachine-precision mass conservation. Given a volume fraction field, reconstructed interface and solenoidal velocityfield, we can solve the equation 15 using an operator splitting algorithm consisting of an x − sweep and a y − sweepfollowing [13]. In order to employ an operator splitting algorithm, the advection of the interface (equation.15), using ∇ · u = 0 , can be written as: ∂C∂t + ∇ · ( u C ) = C ( ∇ · u ) . (16)The above form of the advection equation is essential for performing volume conserving x − direction and y -directionsweeps separately (see [7]). Given a volume fraction ( C ni,j ) and velocity field ( u ni +1 / ,j , v ni,j +1 / ) at the nth time step,the discretised equation 16 is given by, C n +1 i,j = C ni,j + ∆ tr i,j ∆ x (cid:0) δV i − / ,j − δV i +1 / ,j (cid:1) + ∆ t ∆ y (cid:0) δV i,j − / − δV i,j +1 / (cid:1) + C ni,j (cid:18) ∆ tr i,j ∆ x (cid:16) r i +1 / ,j u ni +1 / ,j − r i − / ,j u ni − / ,j (cid:17) + ∆ t ∆ y (cid:16) v ni,j +1 / − v ni,j − / (cid:17)(cid:19) (17)where δV i +1 / ,j = ( ruC ) ni +1 / ,j is the amount of volume fraction fluxed through the right cell face. Similarly, fluxes δV i − / ,j , δV i,j +1 / and δV i,j − / can be computed for other cell faces.Using operator splitting, we can split the above equation as following: C ∗ i,j = C ni,j + ∆ tr i,j ∆ x (cid:0) δV i − / ,j − δV i +1 / ,j (cid:1) + C ∗ i,j (cid:18) ∆ tr i,j ∆ x (cid:16) r i +1 / ,j u ni +1 / ,j − r i − / ,j u ni − / ,j (cid:17)(cid:19) (18) C n +1 i,j = C ∗ i,j + ∆ t ∆ y (cid:0) δV i,j − / − δV i,j +1 / (cid:1) + C ∗ i,j (cid:18) ∆ t ∆ y (cid:16) v ni,j +1 / − v ni,j − / (cid:17)(cid:19) (19)where C ∗ i,j is the intermediate value of the volume fraction. An implicit scheme is used in the first direction and anexplicit scheme in the second direction to maintain the conservation of volume fraction[14]. The order of sweep ofdirection is alternated every timestep [15]("Strang spliting") to achieve second order accuracy in time.The volume flux through cell faces, δV cell − face , is computed geometrically. Consider the schematic in figure 5, wherethe shaded region shows the volume of fluid- in the cell to be fluxed through the right face ( δV i +1 / ,j ). Consideringthe face velocity ( u i +1 / ,j ) to be positive, the flux can be computed as, δV i + ,j = ( ru ) i +1 / ,j V πr ∆ r ∆ y (20)where V is the volume of fluid fluxed through the right face (shown as the shaded region in figure 5), ∆ r is thedistance in the radial direction which contains the volume advected in this timestep and r is the distance to the center ofthis volume from the axis of symmetry. We can calculate ∆ r by considering the conservation of volume fluxed throughthe right face and solving the resulting quadratic equation, which yields, ∆ r = r i + ,j − (cid:113) r i + ,j − r i + ,j u e ∆ t .Using the section of the piece-wise reconstructed interface lying in the volume to be fluxed through the cell-face over ∆ t time-step and employing the Gauss area formula, given by equation 2, we can calculate the volume cut by thisregion.The above small correction in computing ∆ r along with the accurate Gauss formula for axisymmetric simulationsallows us to improve upon the existing volume fraction advection schemes. Existing schemes modify the velocity in D algorithms by using r u for velocity field and use D geometric advection scheme which results in a third ordererror ( O (∆ t h ) , where h is the grid size and ∆ t is the timestep) in mass conservation. We illustrate this by considering x − direction advection of a small volume of fluid through the right face of the cell with a velocity u e , shown as theshaded region in the figure 6 .The existing schemes compute the volume as V c = 2 π ( r e − u e ∆ t )( u e ∆ t )∆ y , where r e is the distance of the eastcell face from the axis of symmetry, ∆ t is the timestep, and ∆ y is the height of the cell. Whereas, the proposed7 PREPRINT - F
EBRUARY
2, 2021Figure 5: The fluxed volume through the right face of the cell when u i +1 / ,j is positive.Figure 6: The shaded region is the volume advected through the east cell face in a single timestep.scheme yields the exact volume, V = 2 π ( r e − ∆ r )∆ r ∆ y with ∆ r = r e − (cid:112) r e − r e u e ∆ t . Thus, the error in volumecalculation is given by, E = π ( u e ∆ t ) ∆ y .We validate the proposed modifications with the following test cases and compare with the results obtained using theopen source multiphase flow solver Gerris[12]. In this test case, a torus of radius . is initialised at (0 . , . in a computational domain of size (2 . , . . Thetorus is advected under the steady state velocity of u = 0 . /r for r > . and v = 0 , where r is the distance from theaxis of symmetry. The fluid is advected timesteps forward in time and then the velocity is reversed to compute timesteps backwards in time. The grid size is / , the grid Courant number(CFL) is chosen to be whichcorresponds to a time step of ∆ t = 0 . . As seen from the figure 7, after timesteps the torus is highlycompressed during the advection as less area (due to axisymmetry) occupies the same volume as we move away fromthe axis of symmetry. We note that the final interface shape matches very well with the initial position of the torus, thusvalidating our algorithm.The relative error in the volume between the initial and final distribution of fluid- for various number of forward andbackward advection time steps is given in table 2. The corresponding relative change in the volume obtained for thesame test case simulated using Gerris flow solver are also given for comparison. We note that the error obtained fromthe present schemes are highly accurate in comparison to those obtained from Gerris.8 PREPRINT - F
EBRUARY
2, 2021 r y At t = 0After 1000 t After 2000 t Figure 7: Interface shape after time steps forward and backward after advecting the torus of radius . with a gridCourant(CFL) number of .Table 3: Results for relative error in the volume for a torus of radius . advected in radial direction forward andbackward in time for different number of timesteps.Relative error in volumeNumber of timesteps Current Solver Gerris Solver . × − . × −
10 2 . × − . × −
100 4 . × − . × − . × − . × − As suggested by Kothe et al.[16], simple linear advection test cases do not reveal the efficacy of advection algorithmsappropriately. Thus, we further test the efficacy of the algorithm by subjecting it to a more severe test case of advectionof a torus in a Hill’s vortex. This is axisymmetric equivalent of the circle in a vortex test case for 2D cases [16]. For thisvelocity field, the interface undergoes strong topological changes including fragmentation and merging due to strongshear effects. Here we use a modified form of Hills’s vortex with a superimposed radial flow field. A torus of radius . is initialised at (0 . , . in a computational domain of size (1 . , . with L = 0 . . The fluid is advected under highlystrained steady state velocity field given by u = 0 . (cid:16) rL ( y − L ) L (cid:17) + . r (21) v = 0 . (cid:20) − (cid:16) y − LL (cid:17) − (cid:0) rL (cid:1) (cid:21) . (22)The fluid is advected timesteps forward in time and then the velocity is reversed to advect timestepsbackwards in time. The grid size is / and the time step is . × − .As seen from figure 8 after timesteps the shape of the interface is highly distorted due to highly strained velocityfield. The final interface shape matches very well with the initial position of the toroid thus validating our algorithm.The relative error in change in volume between the initial and the final distribution of fluid for various number ontime steps is given in table 3. The corresponding relative change in volume for the same test case in Gerris flow solver.We note the error in the proposed scheme, even for larger number of timesteps, is an order smaller compared to theresults from Gerris flow solver. 9 PREPRINT - F
EBRUARY
2, 2021 r y At t = 0After 4000 t After 8000 t Figure 8: Interface shape after time steps forward and backward after advecting the torus placed in a vortex. Therelative error in change in volume is . × − Table 4: Results for relative error in volume for a torus of radius . placed in a complex velocity field advected forwardand backward in time for different number of timesteps.Relative error in volumeNumber of timesteps Current Solver Gerris Solver . × − . × −
10 5 . × − . × −
100 2 . × − . × − . × − . × − . × − . × − In this test case we implement the VOF algorithm presented in this paper for a more complex flow. We solveNavier-Stokes equations in one fluid form given by: ρ ( C ) (cid:18) ∂ u ∂t + u . ∇ u (cid:19) = −∇ p + ∇ · (cid:2) µ ( C ) (cid:0) ∇ u + ( ∇ u ) T (cid:1) (cid:3) + ρ ( C ) g + f γv (23)where u and p are the velocity vector and pressure, respectively, ρ ( C ) and µ ( C ) are the fluid density and viscositywhich are a function of void fraction field, C . We use Chorin’s projection method [17] to solve the above equation 23where we discretise the advection term using a second order ENO scheme [18] and the diffusion terms using centraldifferencing. Surface tension forces, f γv are acting only at the interface and have been modeled as volumetric bodyforce using the continuum surface force model of Brackbill, Kothe, and Zemach [19]. The interface is captured usingCLSVOF algorithm given by Sussman and Puckett [13]. This algorithm is mass conserving and calculates the curvatureand surface normal with high accuracy which is used for surface tension force calculation. The interface is advected bysolving the advection equations for level-set function, φ , and volume fraction, C . We initialize a toroidal bubble ofradius . at (0 . , . in a computational domain of unit size, × . The bottom boundary has an inlet velocity ofunity in the upward axial direction and the right boundary has outflow boundary conditions. The top boundary acts as arigid wall with no-slip and impermeable surface. The density ratio and viscosity ratio is with the Laplace numberof the bubble, La = ρDσµ = 0 . . The incoming axial velocity drags the bubble and flattens it against the top wall,stretching it in the axial direction considerably. Even though the fluid interface undergoes drastic change in its shape,the volume is conserved with a high degree of accuracy with relative error in the volume of . × − .10 PREPRINT - F
EBRUARY
2, 2021 r y At t = 0.0Intermediate time of .2At t = 1.0 Figure 9: Interface shape after time t = 1 . The interface is flattened against the top wall and stretched due to theunderlying velocity. Even when the cross-sectional area changes the volume of the toroidal bubble is maintained veryaccurately. In the present work, we have presented several improvements for the implementation of volume of fluid method inaxisymmetric coordinates. We have presented analytical relations for the reconstruction of piecewise linear interface inaxisymmetric coordinates similar to those given by Scardovelli and Zaleski[10] for cartesian coordinates. The proposedscheme substantially reduces the computational cost in comparison to the iterative schemes usually employed for thereconstruction. Further, we showed that even for axisymmetric coordinate system, machine-precision advection ofvolume fraction field can be achieved. We illustrated the improvements by comparing the results with the popular opensource multiphase flow solver Gerris. Finally, we would like to note that similar modifications in the advection schemefor volume of fluid method in other curvillinear coordinate systems(such as elliptic coordinates) can be derived usingthe approach presented in this work.
References [1] Andrea Prosperetti and Grétar Tryggvason.
Computational methods for multiphase flow . Cambridge universitypress, 2009.[2] Grétar Tryggvason, Ruben Scardovelli, and Stéphane Zaleski.
Direct numerical simulations of gas–liquidmultiphase flows . Cambridge University Press, 2011.[3] Cyril W Hirt and Billy D Nichols. Volume of fluid (vof) method for the dynamics of free boundaries.
Journal ofcomputational physics , 39(1):201–225, 1981.[4] Yuriko Renardy and Michael Renardy. Prost: a parabolic reconstruction of surface tension for the volume-of-fluidmethod.
Journal of computational physics , 183(2):400–421, 2002.[5] David L Youngs. Time-dependent multi-material flow with large fluid distortion.
Numerical methods for fluiddynamics , 1982.[6] James Edward Pilliod.
An analysis of piecewise linear interface reconstruction algorithms for volume-of-fluidmethods . U. of Calif., Davis, 1992.[7] James Edward Pilliod Jr and Elbridge Gerry Puckett. Second-order accurate volume-of-fluid algorithms fortracking material interfaces.
Journal of Computational Physics , 199(2):465–502, 2004.[8] Ruben Scardovelli and Stephane Zaleski. Interface reconstruction with least-square fit and split eulerian–lagrangianadvection.
International Journal for Numerical Methods in Fluids , 41(3):251–274, 2003.[9] David L Youngs. An interface tracking method for a 3d eulerian hydrodynamics code.
Atomic Weapons ResearchEstablishment (AWRE) Technical Report , 44(92):35, 1984.11
PREPRINT - F
EBRUARY
2, 2021[10] Ruben Scardovelli and Stephane Zaleski. Analytical relations connecting linear interfaces and volume fractions inrectangular grids.
Journal of Computational Physics , 164(1):228–237, 2000.[11] Simone Bná, Sandro Manservisi, Ruben Scardovelli, Philip Yecko, and Stéphane Zaleski. Vofi—a library toinitialize the volume fraction scalar field.
Computer Physics Communications , 200:291–299, 2016.[12] Stéphane Popinet. Gerris: a tree-based adaptive solver for the incompressible euler equations in complexgeometries.
Journal of Computational Physics , 190(2):572–600, 2003.[13] Mark Sussman and Elbridge Gerry Puckett. A coupled level set and volume-of-fluid method for computing 3dand axisymmetric incompressible two-phase flows.
Journal of computational physics , 162(2):301–337, 2000.[14] Elbridge Gerry Puckett, Ann S Almgren, John B Bell, Daniel L Marcus, and William J Rider. A high-orderprojection method for tracking fluid interfaces in variable density incompressible flows.
Journal of computationalphysics , 130(2):269–282, 1997.[15] Gilbert Strang. On the construction and comparison of difference schemes.
SIAM journal on numerical analysis ,5(3):506–517, 1968.[16] William Rider and Douglas Kothe. Stretching and tearing interface tracking methods. In , page 1717, 1995.[17] Alexandre Joel Chorin. Numerical solution of the navier-stokes equations.
Mathematics of computation ,22(104):745–762, 1968.[18] Chi-Wang Shu and Stanley Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes.
Journal of computational physics , 77(2):439–471, 1988.[19] Jeremiah U Brackbill, Douglas B Kothe, and Charles Zemach. A continuum method for modeling surface tension.