Front tracking in recirculating flows: a comparison between the TVD and RCM methods in solving the VOF equation

TRANSPORT PHENOMENA AND UNIT OPERATION

Front tracking in recirculating flows: a comparison between the TVD and RCM methods in solving the VOF equation

L. F. Lopes Silva; C. E. Fontes; P. L. C. Lage^*

Programa de Engenharia Química, COPPE, Universidade Federal do Rio de Janeiro (UFRJ), Phone: +(55) (21) 2562-8304, Fax: +(55) (21) 2562-8300, PO Box 68502, ZIP Code 21945-970 Rio de Janeiro - RJ, Brazil. E-mail: [e-mail]

ABSTRACT

The modeling and simulation of mold filling must include a method to capture the interface formed between the inlet fluid and the fluid that was initially in the mold. A commonly used front-capturing method in a Eulerian mesh is the volume-of-fluid (VOF) method. The VOF advection equation solution may show numerical diffusion and/or dispersion and high-order numerical schemes, such as the TVD schemes with dimensional splitting, have to be employed to discretize the convective terms. The present contribution explores the use of RCM for solution of the VOF color-function equation during mold filling with recirculating flows. The Navier-Stokes equations are solved by a finite-volume method using the SIMPLER algorithm. Filling simulations using the TVD and RCM methods are compared. RCM was able to generate diffusion-free results, sharply defining the interface, even when topological changes (generation of droplets) occur.

Keywords: Interface capturing; Mold filling; RCM; TVD schemes; VOF schemes; Computational fluid dynamics.

INTRODUCTION

A variety of industrial processes include a step of filling some kind of cavity (thermoplastic injection molding, reaction injection molding, casting of metal pieces, food extrusion processes, etc). The design of these cavities involves a very expensive trial-and-error procedure in order to avoid filling problems, like heterogeneity in the final piece or air entrapment. This is especially true in the polymer industry. Due to this fact, during the last two decades many researchers have been developing computational methods that allow simulation of cavity filling, which shortens the expensive empirical stage. One of the most challenging problems in constructing a successful simulator is how to capture and track the interface formed between the inlet fluid and the fluid that was initially inside the cavity.

One way of capturing the moving interface between two immiscible fluids is to solve the fluid-flow equations for both fluids in the same Eulerian mesh. The two fluids are modeled as a single continuum with discontinuous properties at the interface. The interface propagation is then obtained through solution of a transient scalar advection equation in the Eulerian mesh. This solution is prone to numerical diffusion and dispersion problems, which are inherent in the numerical solution of hyperbolic equations by any discretization scheme when there are discontinuities in the solution. A fluid-fluid interface is a contact discontinuity, and its resolution demands additional care. In this paper, two numerical solutions to this problem are presented and compared: the use of a high-order TVD scheme (total variation diminishing) developed by Chakravarthy and Osher (1985) and the use of a RCM scheme (random choice method) adapted from Toro (1999). Both schemes use dimensional splitting.

TVD methods are one of the most significant achievements in the development of numerical methods for partial differential equations in the last 20 years or so, although the preliminary ideas can be traced as far back as 1959 to the pioneering work of Godunov, continued later by van Leer. In 1965, Glimm introduced RCM and Colella (1982) contributed to developing and increasing knowledge of the method, considering its advantages and limitations.

In this work, the effectiveness of TVD and RCM for cavity-filling simulation are compared for both nonrecirculating and recirculating flows. In the filling, it is assumed that both fluids had the same properties. Thus, the pressure and velocity fields could be obtained for the steady-state flow throughout any cavity and then used for the filling simulation with both methods. This eliminates any flow field disturbance that would arise in a filling simulation for unsteady fluid flow. These steady-state flow solutions are obtained by solving the transient mass and momentum conservation equations using a finite-volume method and the SIMPLER algorithm until the steady-state flow reached a given tolerance. The effects of surface tension at the interface are not included in the analysis.

FLUID DYNAMICS EQUATIONS

The 2D-Cartesian mass and momentum conservation equations for a Newtonian fluid with constant density and viscosity, used in the simulation of flow through 2D-cavities, are presented below in a dimensionless form.

where

and t is the time; x and y are the spatial coordinates; r is the density; m is the viscosity; n is the kinematic viscosity; u and v are the velocity components in x and y directions, respectively; p is the pressure; g_i is the body force in i direction; D_in is the inlet width and v_med is the mean inlet velocity. The subscript ref refers to the inlet fluid properties.

After applying a finite-volume discretization on a staggered grid using the power-law scheme, the steady-state solution is achieved using the SIMPLER scheme (Patankar, 1980). The numerical code has been tested against benchmark results for the backward facing step flow (Gartling, 1990) with good agreement (Fontes et al., 1999).

The numerical results were integrated using a convergence criterion based on a mixed tolerance that equals the absolute tolerance (atol) added to the relative tolerance (rtol) multiplied by the variable modulus (atol+|x|.rtol<|Dx|). If the absolute tolerance assumes a very small value, the relative tolerance controls the convergence criterion. Absolute and relative tolerances of 10^-6 for the velocity components and of 10^-5 for the pressure were used. The criterion used to assume that the steady state was achieved was the same mixed tolerance criterion described above, but using absolute and relative tolerances of 10^-4 for all the variables in every control volume between two consecutive times.

VOF-TVD SCHEME

The volume-of-fluid (VOF) method is an interface-capturing method (Ferziger & Pèric, 1997). The fluids on both sides of the interface are marked by either massless particles or an indicator (color) function. In the VOF method the advection of the color function F, given by Equation (10) in dimensionless form, has to be solved. The F field values are stored at in the staggered finite-volume grid similarly to the pressure field.

As seen above, U and V are known from the steady-state solution. The simulated physical situation corresponds to the steady flow of a fluid through the cavity that suddenly becomes colored and starts to displace the colorless one. Obviously, if the physical property remains the same, so do the velocity and the pressure fields. For a given volume, F is equal to 0 when the volume has no colored fluid, but it is equal to 1 when the volume is completely filled with colored fluid. The interface has F values between 0 and 1.

Typical problems associated with the solution of hyperbolic equations like Equation (10) (numerical diffusion and dispersion) are reduced by the use of a total variation diminishing (TVD) scheme. Good summaries of TVD schemes are given by Sweby (1984) and Toro (1999). A numerical method is TVD if Equation (11), where x is the transported variable and TV(x) is its total variation, is valid for scalar conservation laws.

Chakravarthy and Osher (1985) described one-dimensional second and third-order accurate TVD schemes with low truncation errors. These schemes were used to solve advective equations, like Equation (10). Goodman and LeVeque (1985) proved that there is no multidimensional TVD scheme, but Chakravarthy and Osher (1985) showed extremely good numerical results using their scheme. Based on its simplicity and good results, the Chakravarthy and Osher (1985) third-order scheme is used in the present work. In order to try to keep the TVD characteristics of the 1-D scheme, a dimensional splitting is applied during the numerical solution of Equation (10). Therefore, the explicit discretized form of Equation (10), with dimensional splitting, is given by Equations (12) and (13), where Dt^1/2 = Dt/2. The splitting algorithm, in accordance with Toro (1999), is described below.

i) For each time step Dt, integrate first in X direction up to Dt^1/2 using Eq. (12).

ii) Then, using F^1/2, integrate in the Y direction up to Dt by Eq. (12), determining F.

iii) In order to achieve second-order accuracy for every-other time step, invert the X and Y integrations in the next time step Dt.

The F fluxes at the faces of the volumes are determined using the upwind-TVD scheme. For example, for the east face, for each j line, the flux is given by Eqs. (14) to (21). In these equations the index P stands for position (i, j), E for position (i+1, j), e for position (i+1/2, j), ee for position (i+3/2, j), w for position (i-1/2, j), n for position (i, j+1/2) and s for position (i, j-1/2). When g equals 1/3 and b equals 4 these equations correspond to a third-order upwind scheme. Equation (14) represents the flux given for the Engquist-Osher first-order scheme (Osher & Chakravarthy, 1984). The minmod function is defined by Equation (19) and it is used to define the limited fluxes given in Equations (17) and (18). The fluxes for the other faces are obtained by changes in subscripts.

Osher and Chakravarthy (1984) prove that their scheme is stable only for the Courant number, given by Eq. (22) for both X and Y directions, equal to or less than 0.4.

Boundary conditions for this equation are as follows: the cells outside the mold inlet are considered to be completely filled with colored fluid (F = 1), all cells outside the mold at the outlet are considered to be completely filled with colorless fluid (F = 0) and there is no F flux at the other boundaries. At the beginning of the filling process, the mold is filled with colorless fluid (F = 0).

VOF-RCM SCHEME

Glimm (1965) introduced the random choice method but it was Colella (1982) who really developed this numerical method, discussing its advantages and limitations. Essentially, RCM is applicable to scalar problems in any number of dimensions without loosing its capacity to maintain the perfect resolution of contact discontinuities with no diffusion or dispersion. As F in Equation (10) is a scalar variable, RCM perfectly matches the necessities of the VOF method. Silva and Lage (2001) show a successful 1-D application of this scheme.

RCM chooses randomly among possible solutions of local Riemann problems associated with the dimensional splitting of Equation (10). The method used in this work is one given by Toro (1999) for a staggered grid in which a local 1-D Riemann problem is solved twice for each direction at each integration step: once between tⁿ and t^n+1/2 and once between t^n+1/2 and tⁿ⁺¹. It is exactly because RCM solves the local Riemann problem that it does not show numerical diffusion. RCM stability is also dependent on the local Courant number, given by Equation (22). For the RCM used in this work, both Courant numbers have to be equal to or less than 1.

The local Riemann problem for Eq. (10) uses the local states (F_L and F_R) around the interface between two control volumes. The exact solution of the local 1-D Riemann problem is shown in Fig. 1-a. It consists in a single contact discontinuity propagation from the origin with velocity u and slope l = 1/u. Naturally, the solution values are equal to F_L on the left of the straight line with slope l = 1/u and equal to F_R on the right.

In order to explain this method, let’s consider its application to the integration of the F value in a given 1-D volume i of a uniform mesh in the X direction, as illustrated in Fig. 1-b, for a time step Dt. First, the local 1-D Riemann problems are solved at both volume interfaces in the first half-time step, Dt^n+1/2 = Dtⁿ/2. At interface (i-1/2), the discontinuity propagates a distance d = X_i-1/2 + C_X,i+1/2DX/2 during Dt^n+1/2. The position r = X_i-1/2 + qⁿDX is sampled by a random number qⁿ Î [0,1]. If r < d, then , otherwise . A similar procedure is applied to the right interface using the same random number qⁿ to obtain . Then, another random number qⁿ⁺¹ is generated, and a similar procedure is applied to obtain the solution for volume i at the end of the time interval, . Thus, two new random numbers are generated for each time step and for each direction. The complete algorithm can be described as follows:

i) Solve the Riemann problems RP and RP to find solutions and ;

ii) For the first half-time step Dt^n+1/2, random sample these solutions based on random number qⁿ: = and =(qⁿ);

iii) Solve Riemann problem RP(,) to find solution ;

iv) For the second half-time step Dt^n+1/2, random sample this solution based on a new random number qⁿ⁺¹: =(qⁿ).

The same dimensional splitting applied in the VOF-TVD scheme is applied in the VOF-RCM scheme, which uses twice the 1-D scheme just described. As RCM is intrinsically a statistical scheme, front behavior may vary between simulation runs for the same problem, which introduces some uncertainty regarding the actual position of the fluid interface. Due to this statistical nature, uncertainty at the fluid interface decreases as the grid is refined. This occurs because more random numbers are generated to obtain the solution at a given time due to the Courant number limitations. The random number generator used was the IMSL DRNUN subroutine that uses the generalized feedback shift register (GFSR) method.

UNIDIMENSIONAL TESTS

In order to verify the TVD and RCM implementations, a one-dimensional test which tracks the movement of a square pulse is used. The simulations were carried out in a uniform grid, using TVD in Fig. 2 and RCM in Fig. 3. Wave positions and shapes are verified after 1.5 s, 3.5 s and 8 s. The agreement with the analytical solution is remarkable if the stability range of the Courant number is respected.

Fig. 2-a shows that the use of coarse grids causes some numerical diffusion in the TVD solution. It was verified that a pulse propagation with a coarse grid deteriorates as the Courant number decreases. When the grid is refined for the same Courant number (Fig. 2-b), this numerical diffusion decreases. It also decreases slightly when the Courant number is increased up to the its maximum (the stability limit) for the same grid, but the results are very similar for C = 0.1 (Fig. 2-a) and 0.4 (Fig. 2-c).

Fig. 3 shows that all RCM results perfectly resolve the contact discontinuity of the interface. However, comparing Figs. 3-a and 3-b it is clear that some inaccuracy exists in the propagated pulse position for coarse grids. This results from errors in the pulse propagation velocity due to the random character of the method. For the finer grid (Fig. 3-b), these errors basically disappear because of the large increase in the quantity of random numbers used to obtain the solution due to the Courant number limitation.

These errors in the position of the contact discontinuities are equivalent to numerical diffusion. This can be seen Fig. 3-c, which shows the mean results of one hundred simulations of the pulse tracking in the coarse grid with different seeds in the random number generator. The behavior of the shape of the propagated pulse is similar to that usually achieved in a single simulation when a diffusive scheme is used. In other words, the uncertainty in the positions of the discontinuities is numerically equivalent to the numerical diffusion of a low-order scheme.

TWO-DIMENSIONAL SIMULATION OF CAVITY FILLING

In order to compare the effectiveness of the TVD and RCM methods associated with the VOF scheme without any perturbation in the flow field, we decided to study the cavity filling using the steady-state flow solution. Two 2-D cavities are used, as shown in Figs. 4-a and 4-b. Dimensions of the mold cavity used were LY = 20 cm and LX = 30.3 cm. The dimensions of the backward facing step (BFS) cavity were 0.1m × 3.0 m.

In a cavity such as the one shown in Fig. 4-a, referred to from now on as mold cavity, a non-recirculating flow with a low Reynolds number Re=0.16 was analyzed. The inlet and outlet widths of the mold are D_in = D_out = 0.1LX. The velocity profile at the inlet is parabolic with a mean velocity of 0.031 m/s (v_med) and the same velocity profile is used at the outlet. Other walls are considered impenetrable and non slip conditions are imposed near all walls. The flow in the BFS cavity (Fig 4-b) with Re=400 has the typical recirculating flow with two opposing vortices. The velocity profile at the inlet is also parabolic, but the mean velocity is equal to 1 m/s. Due to the dimensions of this cavity (LY=30LX), it is possible to apply continuity conditions for the velocities and the pressure at its outlet (Gartling, 1990). The inlet width is D_in = 0.5LX and the outlet width is LX itself. Impenetrable walls and non slip conditions are also imposed on the other walls of this cavity. Simulations use a 200 × 200 grid for the mold cavity and a 200 × 500 grid for the BFS cavity. A time step of 10^-3 s was adequate for all TVD simulations and the BFS-RCM simulations. However, the mold-RCM simulations needed a longer time step (10^-1 s) due to a loss of accuracy of the RCM for low Courant numbers, which is explained in the next section (see the remarks concerning Fig. 9).

As mentioned in Section 2, a power-law scheme is used. The grid Peclet number is defined by Equation (23), where ?the grid length is in the appropriate direction (X or Y).

It is well known that power-law schemes allow numeric diffusion when the grid Peclet number exceeds 10 and the flow is oblique in relation to the grid (Patankar, 1980). In fact, in the BFS simulations, between the recirculating zones, there are Peclet numbers as high as 50 or 60. However, this phenomenon does not invalidate the conclusions of this work because our goal is to compare two surface tracking methods under the same fluid dynamic conditions. Furthermore, as observed in Section 2, velocity and pressure profiles had been successfully tested against Gartling's benchmark.

SIMULATION RESULTS

Figs. 5 and 6 compare the results achieved with the VOF-TVD scheme and the VOF-RCM scheme, respectively, for the nonrecirculating flow inside the mold cavity. The filling patterns are very similar, but a remarkable difference has to be pointed out: as filling continues the TVD scheme loses interface sharpness due to a growing diffusion effect, while the RCM scheme does not show any diffusion. It can be observed that the TVD solution shows the interface deeper inside the cavity than the RCM solution, although the same amount of fluid has entered the mold in both cases. Observing the behavior of the solutions near the bottom wall, it becomes clear that this difference may be explained by the diffusive characteristic of the two-dimensional TVD scheme, which pushes the incoming fluid towards the top wall, allowing the other fluid to remain in the region close to the bottom wall. Due to its statistical characteristic, RCM shows an interface with small fluctuations, which tends to disappear as the filling continues.

Figs. 7 and 8 show the simulation results on the BFS cavity. It is important to observe that the VOF-RCM and VOF-TVD methods show very similar results for the filling of this cavity. However, specifically at the vortex regions, the VOF-RCM method is able to describe the droplets that are pulled away from the colored fluid by the rotating colorless fluid. The VOF-TVD scheme is unable to describe this process. Furthermore, it shows intense numerical diffusion in these regions, which can be seen in Figs. 7-b, 7-c and 7-d.

It is important to point out that Courant number restrictions remain in these two-dimensional simulations. As a matter of fact, Fig. 9-b shows that the VOF-RCM scheme is unable to track the interface when Courant numbers are less than 10^-3. It produces unrealistic results. On the other hand, the VOF-TVD scheme (Fig. 9-a) does not show the same problem and tracks the interface perfectly.

It is important to point out that all the simulations described above show an error in mass conservation of less than 5%, except for the simulations using the VOF-RCM scheme shown in Fig. 9-b.

Although the VOF-RCM scheme resolves the interface perfectly, it must be pointed out that with different seeds in the random number generator, different simulated results are obtained, even though their overall filling patterns are similar. Thus, both VOF-TVD and VOF-RCM schemes seem to be adequate to describe the overall picture of the filling pattern in the mold and the BFS cavities.

CONCLUSIONS

It was shown that the VOF-RCM scheme is able to generate diffusion-free results, sharply defining the interface, even when topological changes (droplet generation) occur due to intense recirculating flow. However, the statistical nature of RCM demands fine grids in order to simulate the interface position correctly. On the other hand, the VOF-TVD scheme shows some numerical diffusion but is less sensitive to low Courant numbers than the VOF-RCM scheme. It is important to recognize that the front positions indicated by the RCM method may vary from one simulation to the next due to the statistical nature of the method. It was shown that these uncertainties are equivalent to numerical diffusion.

ACKNOWLEDGEMENTS

The authors acknowledge the financial support obtained from CNPq under grant number 520660/98-6.

NOMENCLATURE

C_i Courant number in direction I (-) D_in inlet width, M D_out outlet width, M F color function (-) Riemann problem solution (-) g_i body force in i direction, m/s² P pressure, Pa P adimensional pressure (-) Re_ref inlet Reynolds number (-) t time, s T stress tensor, Kg/m^.s² U velocity component in x direction, m/s V velocity component in y direction, m/s v_med mean inlet velocity, m/s Y spatial coordinate, m X spatial coordinate, m b, g coefficients in upwind scheme (-) qⁿ random number Î [0,1] (-) G adimensional viscosity (-) j adimensional density (-) df_i limited flux for face I (-) r fluid density, Kg/m³ r_ref inlet density, Kg/m³ m fluid viscosity, Pa^.s m_ref inlet viscosity, Pa^.s t adimensional time (-) x transported variable (-) D deviation associated with a quantity (-)

Received: February 8, 2003

Accepted: August 12, 2004

• Chakravarthy, S.R. and Osher, S., A New Class of High Accuracy TVD Schemes for Hyperbolic Conservation Laws, AIAA paper 85-0363 (1985).
• Colella, P., Glimms Method for Gas Dynamics, SIAM J. Sci. Stat. Comput., vol. 3, no. 1, pp. 76-110 (1982).
• Ferziger, J.H and Pèric, M., Computational Methods for Fluid Dynamics, Springer, New York (1997).
• Fontes, C.E., Lage, P.L.C. and Pinto, J.C., Filling of Rectangular Molds Using the VOF Method and a TVD Discretization, Proceedings of the XV Brazilian Congress of Mechanical Engineering (Cobem), Águas de Lindóia, São Paulo, Brazil (1999).
• Gartling, D.K., A Test Problem for Outflow Boundary Conditions - Flow over a Backward Facing Step, International Journal for Numerical Methods in Fluids, vol. 11, pp. 953-967 (1990).
• Glimm, J., Solution in the Large for Nonlinear Hyperbolic Systems of Equations. Comm. Pure. Appl. Math., vol. 18, pp. 697-715 (1965).
• Godunov, S.K., A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics, Mat. Sb., vol. 47, pp. 357-393 (1959).
• Goodman, J.B. and LeVeque, R.J., On the Accuracy of Stable Schemes for 2-D Scalar Conservation Laws, Math. Comp., vol. 45, p. 15 (1985).
• Osher, S. and Chakravarthy, S., Very High Order Accurate TVD Schemes, ICASE Report no. 84-44, ICASE NASA Langley Research Center, Hampton, Virginia (1984).
• Patankar, S.V., Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York (1980).
• Silva, L.F.L.R. and Lage, P.L.C., Método de Escolha Aleatória em Malhas Não Uniformes: Solução da Equação da Advecção Unidimensional, 4ş Congresso Brasileiro de Engenharia Química em Iniciação Científica COBEC-IC, Maringá PR (CD-ROM) (2001).
• Sweby, P.K., High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws, SIAM Journal of Numerical Analysis, vol. 21, no. 5, pp. 995-1011 (1984).
• Toro, E.F., Riemann Solvers and Numerical Methods for Fluid Dynamics - A Practical Introduction, 2nd Edition, Springer (1999).
• Van Leer, B., Toward the ultimate conservative difference scheme. V. A second-order sequel to Godunovs Method, Journal of Comp. Physics, vol. 32, pp. 101-136 (1979).

*

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  14 Mar 2005
• Date of issue
  Mar 2005

History

• Received
  08 Feb 2003
• Accepted
  12 Aug 2004