An implicit technique for solving 3D low Reynolds number moving free surface flows

https://doi.org/10.1016/j.jcp.2008.04.017Get rights and content

Abstract

This paper describes the development of an implicit finite difference method for solving transient three-dimensional incompressible free surface flows. To reduce the CPU time of explicit low-Reynolds number calculations, we have combined a projection method with an implicit technique for treating the pressure on the free surface. The projection method is employed to uncouple the velocity and the pressure fields, allowing each variable to be solved separately. We employ the normal stress condition on the free surface to derive an implicit technique for calculating the pressure at the free surface. Numerical results demonstrate that this modification is essential for the construction of methods that are more stable than those provided by discretizing the free surface explicitly. In addition, we show that the proposed method can be applied to viscoelastic fluids. Numerical results include the simulation of jet buckling and extrudate swell for Reynolds numbers in the range [0.01, 0.5].

Introduction

A very active and important research area in Computational Fluid Dynamics is the development of computational methods for free surface flows. This can be explained by the fact that flows with free surfaces are more difficult to resolve in general than confined flows. The position of the boundary is known only at the initial time, and its new position has to be determined as part of the solution. The MAC method (Marker-and-Cell), introduced by Harlow and Welch [18], was one of the first successful attempts to simulate viscous, incompressible, transient flows with free surfaces. MAC was derived from the discretization of the Navier–Stokes equations in primitive variables by finite differences on a uniform staggered mesh. In this method, the shape of the free surface is determined by cells that are partially filled. Variants of the MAC method with distinct versions are widely available in the literature (see e.g. [2], [8], [9], [18], [19], [25], [33], [41], and many others). However, one feature common to almost all these techniques is the explicit time discretization of the momentum equations by Euler’s method. Indeed, for inertial flows (Re>1) explicit methods do not impose too severe a restriction on the time step so that numerical solutions of free surface flows can be obtained in reasonable time (e.g. [38]). However, Reynolds numbers of order 10-1 to 10-3 can be easily found in applications involving the flow of polymers such as extrudate swell, injection moulding and jet buckling. For low Reynolds number flows, explicit methods have a stringent stability restriction leading to very small time steps. One way to overcome this restriction is to use an implicit discretization that, one might expect, would lead to a substantial reduction in CPU time. There would appear to be only a few papers (e.g. [13], [15], [17], [31], [33], [47]) presenting implicit formulations for free surface flows, and these tend to be for two-dimensional flows where the Reynolds number is not excessively small. However, there are some works that employ the level set method to model free surface flows. For instance, Fedkiw et al. [12] presented a numerical method for treating interfaces using an Eulerian scheme for multi-material flows while Nguyen et al. [26] solved the Euler equations to simulate multiphase flows of inviscid fluids. To calculate the velocity field they used the projection method in a similar fashion as in the MAC method [18]. They reported some results for 1D, 2D and 3D problems. One interesting work that deals with free surface viscous flows was presented by Vincent and Caltagirone [45] (see also [46]). They solved the full Navier–Stokes equations for interfacial flows on a staggered grid using the finite volume method. The interface was represented by an advection equation which was treated by a high order TVD (Total Variation Diminishing) scheme. Their results include the simulation of 3D advection of a sphere and the simulation of a 2D Newtonian jet filling a square box.

More recently, the authors presented a novel two-dimensional implicit method (see Oishi et al. [28]) for unsteady two-dimensional free surface flows using a Marker-and-Cell approach. In that method the implicit Euler scheme was employed in the discretization of the diffusion terms while the free surface boundary conditions were discretized implicitly. Numerical experiments obtained by Oishi et al. [28] suggested that when the equation for the pressure on the free surface was discretized explicitly, a parabolic-type stability condition was required to be imposed on the time step, independently of whether the diffusion terms of the Navier–Stokes equations were discretized explicitly or implicitly. Indeed, it was found numerically that in order to obtain a stable implicit solver for the Navier–Stokes equations it was crucial that the pressure condition at the free surface was required to be solved implicitly. However, when the pressure equation at the free surface was discretized implicitly it coupled the velocity and pressure fields so that a much larger linear system had to be solved at each time step. To overcome this, Oishi et al. [28] developed numerical techniques for uncoupling the velocity and pressure fields while still maintaining the implicit discretization of the pressure equation at the free surface. More recently, Oishi et al. [29] have performed a rigorous stability analysis on discretizations of a paradigm problem, the 1-D heat equation, solved on a staggered grid with the implicit boundary conditions replaced by explicit boundary conditions and this model problem has shed light as to why the approach described above is effective.

In this work we employ the ideas presented by Oishi et al. (see [28]) to develop a stable 3D implicit method for solving incompressible free surface flows. This novel formulation combines an accurate projection method (described in [30]) with a new formula for the pressure-update and an implicit technique for computing the pressure on the free surface. This last point is very important because, as was shown in [28], the implicit discretization of the boundary conditions at the free surface is crucial for achieving a stable implicit method. The performance of the numerical method is demonstrated by simulating three-dimensional free surface flow problems with low Reynolds numbers and moving free surfaces. Moreover, we show that the proposed method can also be applied to viscoelastic fluids.

This paper is organized as follows. First the governing equations together with the boundary conditions are described. Section 3 presents the details of the numerical method: a description of the projection method, the implicit formulation for the pressure at the free surface, the algorithm itself, the detailed finite difference approximations and a discussion on the time step stability restriction. Numerical results are given in Section 4. The flow of an Oldroyd-B fluid inside a 3D tube is simulated and compared with the corresponding analytic solution. A comparison between the explicit and implicit techniques is then effected. Finally, the implicit technique is employed to solve complex free surface flows: both extrudate swell of a 3D Oldroyd-B fluid and the buckling of a Newtonian jet are simulated; extremely small Reynolds numbers were used in both cases. Section 5 contains some concluding remarks.

Section snippets

Governing equations and boundary conditions

In dimensionless conservative form, the equations for incompressible viscous flow in primitive variables can be written asut+·(uu)=-p+α1Re2u+β·S+1Fr2g,and the mass conservation equation as·u=0,where S is the non-Newtonian part of the extra stress tensor (hereafter simply called the non-Newtonian stress tensor) that is defined by an appropriate constitutive relationship. For Newtonian flows S=0 so this can be achieved by setting α=1 and β=0 in the momentum Eq. (1). In this work we shall

Numerical method

The numerical method for solving Eqs. (1), (2) will be based on the projection method (pioneered by Chorin [9]) while the constitutive Eq. (3) is approximated by a second-order finite difference method. In summary, the solution of the momentum Eq. (1) is obtained by calculating a provisional velocity field followed by the solution of an elliptic equation to enforce mass conservation (2).

The projection methods are based on the Helmholtz decomposition (see [11]) which states that every smooth

Numerical simulation of low Re free surface flows

The implicit technique described in the previous sections was implemented into the Freeflow3D code (see Castelo et al. [7]) to simulate unsteady three-dimensional low Reynolds number free surface flows of Newtonian and Oldroyd-B fluids. It was used to simulate the flow of an Oldroyd-B fluid in a tube and the numerical results were compared with the analytic solution. Time-dependent extrudate swell and jet buckling for very small Reynolds numbers were also simulated. The results presented in

Concluding remarks

This paper has dealt with the development of an implicit technique for solving three-dimensional free surface flows of Newtonian and Oldroyd-B fluids. The momentum equations were solved by the Crank–Nicolson technique and an implicit method for treating the pressure on the free surface was developed so that the pressure and the velocity fields could be decoupled. The Oldroyd-B constitutive equation has been solved using an explicit finite difference technique developed by Tomé et al. [44]. The

Acknowledgements

The authors would like to acknowledge the financial support of FAPESP (Fundação de Amparo a pesquisa do Estado de São Paulo) (Grants No. 03/12612-9 and 04/16064-9), CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) (Grant No. 300106/2005-0) and CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nı´vel Superior) (Grant No. 136/05-CAPES/GRICES). The second author has been supported by CNPq under grants 474040/2003-8, 523141/94 and CAPES pos-doc Grant No. 0121/07-0.

This work

References (47)

  • H. Miyata et al.

    Finite-difference simulation of non-linear waves generated by ships of arbitrary three-dimensional configuration

    J. Comput. Phys.

    (1985)
  • H. Miyata

    Finite-difference simulation of breaking waves

    J. Comput. Phys.

    (1986)
  • D.Q. Nguyen et al.

    A boundary condition capturing method for incompressible flame discontinuities

    J. Comput. Phys.

    (2001)
  • B. Perot et al.

    A moving unstructured staggered mesh method for the simulation of incompressible free surface flows

    J. Comput. Phys.

    (2003)
  • W.E. Pracht

    A numerical method for calculating transient creep flows

    J. Comput. Phys.

    (1971)
  • F.S. de Sousa et al.

    A front-tracking/front-capturing method for the simulation of 3D multi-fluid flows with free surfaces

    J. Comput. Phys.

    (2004)
  • M.F. Tomé et al.

    GENSMAC: a computational marker-and-cell method for free surface flows in general domains

    J. Comput. Phys.

    (1994)
  • M.F. Tomé et al.

    A numerical method for solving three-dimensional generalized Newtonian free surface flows

    J. Non-Newtonian Fluid Mech.

    (2004)
  • M.F. Tomé et al.

    Die-swell, splashing drop and a numerical technique for solving the Oldroyd-B model for axisymmetric free surface flows

    J. Non-Newtonian Fluid Mech.

    (2007)
  • S. Vincent et al.

    A one-cell local multigrid method for solving unsteady incompressible multiphase flows

    J. Comput. Phys.

    (2000)
  • B. Yang et al.

    A second-order boundary-fitted projection method for free surface flow computations

    J. Comput. Phys.

    (2006)
  • M. Alves et al.

    A convergent and universally bounded interpolation scheme for the treatment of advection

    Int. J. Numer. Methods Fluids

    (2003)
  • A.A. Amsden et al.

    A simplified MAC technique for incompressible fluid flow calculations

    J. Comput. Phys.

    (1970)
  • Cited by (27)

    • A least-squares particle model with other techniques for 2D viscoelastic fluid/free surface flow

      2020, Journal of Computational Physics
      Citation Excerpt :

      The adopted mathematical model is constructed by constitutive differential equations which are usually obtained by continuum approach, and the corresponding micro properties are obtained empirically (see [1]), for example the Oldroyd-B [3–5,8], Giesekus [2] and FENE-P [11] models and so on. Moreover, the shear thinning, overshooting and the elastic instability characters (see [1–15]) of viscoelatic fluid flow are known as the international academic hotspots in recent years. The above rheological behaviors usually appear in the high Weissenberg number problem (HWNP) or viscoelastic free surface flow.

    View all citing articles on Scopus
    View full text