Elsevier

Journal of Computational Physics

Volume 291, 15 June 2015, Pages 151-176
Journal of Computational Physics

hp-Adaptive time integration based on the BDF for viscous flows

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

Abstract

This paper presents a procedure based on the Backward Differentiation Formulas of order 1 to 5 to obtain efficient time integration of the incompressible Navier–Stokes equations. The adaptive algorithm performs both stepsize and order selections to control respectively the solution accuracy and the computational efficiency of the time integration process. The stepsize selection (h-adaptivity) is based on a local error estimate and an error controller to guarantee that the numerical solution accuracy is within a user prescribed tolerance. The order selection (p-adaptivity) relies on the idea that low-accuracy solutions can be computed efficiently by low order time integrators while accurate solutions require high order time integrators to keep computational time low. The selection is based on a stability test that detects growing numerical noise and deems a method of order p stable if there is no method of lower order that delivers the same solution accuracy for a larger stepsize. Hence, it guarantees both that (1) the used method of integration operates inside of its stability region and (2) the time integration procedure is computationally efficient. The proposed time integration procedure also features a time-step rejection and quarantine mechanisms, a modified Newton method with a predictor and dense output techniques to compute solution at off-step points.

Introduction

Time integration of mathematical models describing unsteady phenomena is a key issue in computational sciences and engineering. Many of these models can be cast in the form of an initial value problem in time and a boundary value problem in space. There is a rich literature dealing with methods and algorithms to solve Ordinary Differential Equations (ODE) in time for both non-stiff problems (see Ref. [1]) and stiff problems (see Ref. [2]). For a given ODE, the choice of one method over the others is determined by considerations of stability, accuracy and computational efficiency [3]. Note that a method is said to be more computationally efficient than another if it can achieve the same accuracy for a much lower cost.

This paper focuses on unsteady incompressible viscous flow simulations based on the Navier–Stokes equations. Once space discretization is performed, the resulting time-dependent discrete problem is not a system of ODEs but rather a system of Differential Algebraic Equations (DAE) of index-2. Furthermore, they constitute a set of stiff differential equations. Because of these two properties, time integration of the Navier–Stokes equations suffers order and stability restrictions which greatly limit the range of possible time-stepping schemes that can be used [2].

The most popular methods to perform time integration of the incompressible Navier–Stokes equations are the implicit θ-schemes (especially the backward Euler and Crank–Nicolson (CN) method), the generalized-α schemes, the Gear scheme (or second-order Backward Euler), Runge–Kutta methods and Rosenbrock methods (see for instance Refs. [4], [5], [6], [7], [8]). Note that we only discuss implicit methods because they perform tremendously better than explicit ones for stiff problems. The CN and Gear schemes have been the most used methods in the engineering literature because they are cheap (both in terms of memory requirement and computational time to perform a single step) and their implementation in a computer code is relatively easy. However, as the generalized-α methods, they are at best only second-order accurate so that the computational efficiency is limited to relatively low-accuracy solutions. Indeed, high-accuracy solutions can be obtained at a much lower computational cost by using higher-order methods (see e.g. Refs. [6], [7], [8]). Implicit Runge–Kutta (IRK) methods and Rosenbrock methods can deliver higher convergence rates but at the price of increasing the size of the linear algebraic systems when increasing the order of convergence. Diagonally Implicit Runge–Kutta (DIRK) and Rosenbrock methods suffer significant order reduction and stability limitation when applied to index-2 DAE (see for instance Ref. [7]). While very attractive, Singly Diagonally Implicit Runge–Kutta (SDIRK) methods must be rejected because they cannot achieve orders of convergence higher than those of the CN scheme (i.e. first-order accuracy for the pressure and second-order accuracy for the velocity). Among the IRK methods, the Radau IIA time-stepping schemes have the best stability properties given that they are L-stable [2] for index-2 DAEs so that they are unconditionally stable and they have excellent dissipation properties for unresolved modes in the solution. (Also they are energy-conserving methods [9].) Their only limitation for index-2 DAEs is that they deliver a lower time accuracy for the pressure than for the velocity [10]. Other IRK methods are not L-stable (most are even only conditionally stable) or require a higher number of stages to reach the same order of convergence so that they are less computationally attractive. As an example, the 3-stage Radau IIA is 5th order accurate for the velocity, 3rd order accurate for the pressure and has no stability restriction on the stepsize. However, the method leads to linear systems three times as large as that of the CN or Gear schemes. Thus, their use for two-dimensional calculations is unpractical and is unreasonable for three-dimensional simulations given their memory requirement [8].

All previously mentioned methods are one-step methods with the exception of the Gear scheme. The main advantage of one-step methods is that they are self-starting. However, the Gear time-stepping procedure has been very popular in the engineering literature because it is L-stable and second-order accurate for both velocity and pressure. The Gear method belongs to a class of implicit multistep time integrators known as the Backward Differentiation Formulas (BDF) [11]. An accuracy to cost comparison between the 2- and 3-stage Radau IIA IRK and the BDF of order 1 to 6 (referred to as BDF-1 to BDF-6 respectively in what follows) is provided in Ref. [8]. It has been shown that, when used within their stability region, the BDF methods are much more computationally efficient than the IRK for three-dimensional applications. BDF-1 is the first-order implicit backward Euler method while BDF-2 is the Gear time-stepping scheme. We note that higher-order BDF methods are much less popular because their stability properties deteriorate as the order increases so that BDF methods of orders higher than 6 are unconditionally unstable [1]. Stability properties become even worse for variable stepsize formulas. However, BDF methods are easy to implement and cheap: their cost is almost independent of the order of accuracy since the only additional requirement for higher order methods is to store additional solution vectors at previous instants.

In Ref. [8], it has been demonstrated in the context of Arbitrary Eulerian–Lagrangian formulation that the BDFs can be used to obtain higher order solutions in a cost effective manner. The main restriction arising from the stability limitation of the BDF schemes is that high-order methods (orders 3 to 6) cannot be used to obtain relatively low accuracy solution. This is not a severe restriction because low accuracy solutions can be obtained in a cost effective manner by the low-order BDF methods. Indeed, higher-order BDFs become advantageous only when high accuracy solutions are sought. This means that the stability limitation of BDF methods requires to carefully select which method within the family should be used based on the level of requested accuracy. However, the selection is not easily performed a priori unless it is based on previous experience. It is thus natural to look for an automatic order selection process.

In the field of computational ODEs, the concept of adaptive selection of order and stepsize is not a new idea. Indeed, adaptivity is a powerful approach in computational methods to yield cost effective algorithms. In the literature, adaptive procedures have been devised for almost all time-stepping procedures (see for instance Refs. [12], [13], [14]). These algorithms perform stepsize (h-adaptivity) and/or order selection (p-adaptivity). From them, a set of efficient adaptive solvers have been written. A review of these along with their performances on a set of test cases is presented in Ref. [1], Section III.7 p. 421.

ODE solvers based on multistep methods have been introduced and popularized by the work of Gear [15]. The well-known DIFSUB ODE solver [16] based on Adams methods for non-stiff problems and BDF for stiff problems has been very successful in the 70s.

Later, hp-adaptivity has been used to design efficient algorithms to solve stiff problems based on the BDF (see e.g. Ref. [17]). Among them, DASSL developed by Petzold has proven extremely efficient to solve ODEs and DAEs of index 0 and 1 (see Ref. [18] for a detailed presentation). In this paper we will revisit these techniques to design a cost-effective algorithm for adaptive time integration of the Navier–Stokes equations. The basic ingredients are:

  • A stability indicator that detects growing numerical noise in approximate solutions. Our approach is based on the idea of Shampine and Gordon [19] which requires that higher order terms in the Taylor series expansion of the local truncation error to form a sequence of decreasing magnitude. This allows to adaptively perform order selection based on the stability state of the current numerical solution sought.

  • A local error estimator based on the first term in the local truncation error and an error controller. This allows to adaptively perform on-the-fly stepsize selection based on a user-supplied error tolerance.

  • A mechanism to reject and restart time steps for which the preset error tolerance has not been achieved.

  • A predictor based on extrapolation from the BDF to yield a good initial guess for the Newton iterative process at each time step.

  • A modified Newton method whose convergence criteria (to control the iterative error) are based on the chosen error tolerance (to control the discretization error).

To summarize h-adaptivity controls the solution accuracy while p-adaptivity controls the computational efficiency of the time integration process. The automatic control of the error frees the user from the tedious task of determining and prescribing beforehand time stepsize distribution as a function of time that will yield solution of the prescribed accuracy. The control of the integration order allows to obtain solution in a cost effective manner.

The proposed algorithm is based on ideas similar to those used in previous adaptive BDF codes and especially in DASSL. However, the algorithm presented in this paper differs from previous ones in a number of ways, in particular:

  • The hp-algorithm has been specifically tuned for incompressible flow problems for which the stability restriction of the BDF depends on the Reynolds number (see Section 2.4).

  • We use the variable stepsize formulation of the BDF to achieve a better compromise between hp-adaptivity and the modified Newton method (see Sections 2.2 and 3.4).

  • The stability test is based on the following stability definition for the BDF: a method is stable if there is no lower order methods that deliver the same error level for a larger stepsize (see Section 3.2).

  • A complete mechanism to determine if a time step is to be rejected which also includes a recovery procedure (see Section 3.5).

  • The dense output mechanism based on 6th order accurate spline reconstruction to deliver continuous solution representations over the entire integration interval (see Section 3.6).

The paper is organized as follows. Section 2 summarizes the main properties of the BDF time integrators for ODEs and DAEs. The predictor and corrector formulas are presented; error estimators are derived in Section 2.3 and the stability test introduced in Section 2.4. Section 3 provides a thorough description of the adaptive algorithm. In particular stepsize selection is presented in Section 3.1, order selection in Section 3.2 and the rejecting-step mechanism in Section 3.5. We also discuss the starting procedure, the error controller, the modified Newton method and the dense output technique. Finally, Section 4 gives detailed numerical examples of typical incompressible flow problems. The behavior and performances of the proposed time integration algorithm are illustrated and discussed.

Section snippets

The Backward Differentiation Formulas

In this section, we summarize a number of theoretical elements that are important for the development of the variable stepsize implicit Backward Differentiation Formulas (BDF). These elements are presented in details in a number of books on computational methods for initial-value problems. The interested reader is referred to the following references [18], [3], [1], [2]. In most references the theory is presented for constant stepsize methods while we focus here on variable stepsize methods and

Adaptive procedure

Algorithm 1 provides an overall sketch of the hp-adaptive time-integration procedure. At each time step we determine the order of the BDF based on stability consideration (see Section 3.2). Next, the error controller of the stepsize selection process predicts the largest stepsize to be used based on the accuracy requirement (see Section 3.1). Once the order and stepsize have been determined, the corrector coefficients for the current step are calculated (see Appendix B) and the predictor

Flow over a circular cylinder

We consider the two-dimensional flow over a circular cylinder characterized by the Reynolds number (see Section 2.4). We solve the non-dimensional form of the Navier–Stokes equations using the cylinder diameter D as a reference length and the free-stream velocity U as a reference velocity (so that D/U gives a reference time). Therefore, the Reynolds number is defined as Re=UD/ν. As classically done, we use the dynamic pressure applied on the unit cylinder 12ρU2D as a reference force. Note

Conclusion

This paper has presented an hp-adaptive procedure based on the Backward Differentiation Formulas to yield efficient time integration of the incompressible Navier–Stokes equations. The stepsize selection (h-adaptivity) is based on a local error estimate and an error controller to guarantee that the numerical solution accuracy is within a prescribed user tolerance. The proposed algorithm includes a mechanism to reject and restart steps for which the preset error tolerance has not been met.

The

References (31)

  • U.M. Asher et al.

    Computer Methods for Ordinary Differential Equations and Differential–Algebraic Equations

    (1998)
  • P.M. Gresho et al.

    Incompressible Flow and the Finite-Element Method: Volume 1, Advection–Diffusion

    (1998)
  • A. Montlaur et al.

    High-order time integration for unsteady incompressible flows

    Int. J. Numer. Methods Fluids

    (2012)
  • C.W. Gear

    Numerical Initial Value Problems in Ordinary Differential Equations

    (1971)
  • A. Malidi et al.

    A study of time integration schemes for the numerical modelling of free surface flows

    Int. J. Numer. Methods Fluids

    (2005)
  • Cited by (0)

    View full text