BE-BDF2 Time Integration Scheme Equipped with Richardson Extrapolation for Unsteady Compressible Flows

Abstract

In this work we investigate the effectiveness of the Backward Euler-Backward Differentiation Formula (BE-BDF2) in solving unsteady compressible inviscid and viscous flows. Furthermore, to improve its accuracy and its order of convergence, we have equipped this time integration method with the Richardson Extrapolation (RE) technique. The BE-BDF2 scheme is a second-order accurate, A-stable, L-stable and self-starting scheme. It has two stages: the first one is the simple Backward Euler (BE) and the second one is a second-order Backward Differentiation Formula (BDF2) that uses an intermediate and a past solution. The RE is a very simple and powerful technique that can be used to increase the order of accuracy of any approximation process by eliminating the lowest order error term(s) from its asymptotic error expansion. The spatial approximation of the governing Navier–Stokes equations is performed with a high-order accurate discontinuous Galerkin (dG) method. The presented numerical results for canonical test cases, i.e., the isentropic convecting vortex and the unsteady vortex shedding behind a circular cylinder, aim to assess the performance of the BE-BDF2 scheme, in its standard version and equipped with RE, by comparing it with the ones obtained by using more classical methods, like the BDF2, the second-order accurate Crank–Nicolson (CN2) and the explicit third-order accurate Strong Stability Preserving Runge–Kutta scheme (SSP-RK3).

1. Introduction

To understand a physical phenomenon the mathematical model that mimics it needs often to be discretized since its analytical solution is not available. Therefore, the numerical techniques used to obtain an approximate solution of mathematical models have become an indispensable part of modern science.

When modelling the behaviour of an unsteady flow, the Partial Differential Equations (PDEs) that describe it often lead to a stiff problem, since the physical phenomenon is characterised by a wide range of spatial and temporal scales, e.g., turbulent flows. In this case, implicit methods are potentially a better choice with respect to explicit ones, due to their better stability properties. Nevertheless, although very promising implicit methods have been developed by mathematicians in the past and are well known in the literature, see for example [1,2,3,4] just to name a few, their use within Computational Fluid Dynamics (CFD) codes is somewhat limited. This is probably due to the greater computational complexity of implicit methods with respect to explicit ones and to the greater CPU time and memory requirements needed by each time-step to advance the solution in time, which is a really important aspect when one has to deal with problems that involve a large number of degrees of freedom (DOFs). Nevertheless, starting from implicit “classical” methods, like the Backward Differentiation Formulae (BDF) [5] or the second-order accurate Crank–Nicolson scheme (CN2) [6], many other variants have been developed in the past with the aim to improve their performance. Some examples of methods that have been developed to improve the stability region of the BDF schemes are the MEBDF [2,7,8] or the TIAS [3,9] schemes, that are in fact multi-stage methods A-stable up to order 4 and 6, respectively, or the second derivative method proposed by Enright [1,10] that has obtained a single-stage method A-stable up to order 4 by using a higher derivative of the solution. Another way to improve the stability properties of the BDF schemes was found with the Composite-Backward Differentiation Formulae (C-BDF) [4,11,12], which are schemes that inherit the L-stability property of the second-order accurate BDF scheme and that, to the knowledge of the author, have been mainly used until now for electromagnetic transient simulations [13] and for thermal radiative diffusion problems [14], and only recently has their potential in solving structural mechanics [15] and fluid dynamics [16] problems been investigated. Furthermore, note that in the above cited references, the C-BDF schemes have been used without the Richardson Extrapolation (RE) technique.

In view of this, the purpose of this article is to contribute to minimise the gap between the mathematical formulation and the study of theoretical aspects of the C-BDF schemes and their practical use to approximate unsteady compressible inviscid and viscous flows, even when they are equipped with the RE technique. Please note that, to the best of the author’s knowledge, the coupling of the C-BDF schemes with the RE technique is investigated here for the first time. Among all the possible C-BDF schemes, in this work we consider the implementation in a CFD code, based on a discontinuous Galerkin (dG) method for the spatial approximation of the Backward Euler-Backward Differentiation Formula (BE-BDF2) [11,17] (note that in these references, this scheme was named with the more general acronym C-BDF2). The dG method is a high-order accurate method [18,19,20], which is very well suited to cope with problems requiring a high accuracy solution and, thanks to its compact form, for parallel computations. Furthermore, thanks to its favourable dissipation and dispersion properties, the dG method has been proved to be very effective to solve turbulent flows with both the Direct Numerical Simulation (DNS) [21,22] and the Large Eddy Simulation (LES) [23,24]. The BE-BDF2 scheme is very similar to the TR-BDF2 scheme [4] and can be considered a variant of it. According to the Second Dahlquist barrier [25], this scheme is only second-order accurate, but, as highlighted in this work, this relatively low order of accuracy can be improved thanks to a very simple and efficient numerical technique, known in the literature as Richardson Extrapolation (RE) [26,27]. The RE technique, named by Richardson himself as a “deferred approach to the limit”, is able to dramatically improve the accuracy and raise the order of convergence of any numerical method with a very simple idea. This simple idea is to compute two solutions with the same numerical method, one on a coarse grid and one on a fine grid, and to use a simple linear combination of these two approximations with the goal of eliminating the lowest-order error term(s) of their asymptotic error expansion. The final result of the RE technique is a better approximation of the final solution and an increasing of the order of accuracy of the numerical method used.

To clearly determine and evaluate the effectiveness of the BE-BDF2 scheme, in its standard form and equipped with the RE technique, in this work we present the numerical results of the following test cases: (i) the inviscid isentropic vortex convection; (ii) the laminar vortex shedding behind a circular cylinder. For the first test case, the performance of the numerical schemes presented in this work will be compared with the ones obtained with classical BDF2 and CN2 methods, highlighting that, for a fixed time-step size, the use of these schemes allows us to achieve more accurate solutions, especially when equipped with the RE technique. For the second test case, that is a stiff-problem, the BE-BDF2 and BE-BDF2 & RE schemes will be compared with the explicit SSP-RK3, showing that the force coefficients predicted by all the schemes are very similar, despite the very large time-step used for the implicit schemes, which is about 220 times larger than the one employed by the explicit scheme for stability reason.

In the next Section are reported the Navier–Stokes governing equations. Section 3 briefly describes the main building blocks of the dG spatial discretization. Section 4 is devoted to the description of the BE-BDF2 and RE methods and to the main computational aspects related to their implementation. The numerical results in Section 5 reveal the performance of the BE-BDF2 and BE-BDF2 & RE schemes. Findings, main conclusions and possible future works are summarised in the last Section.

2. The Governing Equations

The equations governing the behaviour of viscous compressible flows, i.e., the Navier–Stokes equations, can be written as:

$$\begin{array}{cc}\hfill & \frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_j\right)=0,\hfill \\ \hfill & \frac{\partial }{\partial t}\left(\rho u_i\right)+\frac{\partial }{\partial x_j}\left(\rho u_i{u}_j\right)=-\frac{\partial p}{\partial x_j}+\frac{\partial {\tau }_{ij}}{\partial x_j},\mathrm{for}i,j=1,..,d\hfill \\ \hfill & \frac{\partial }{\partial t}\left(\rho E\right)+\frac{\partial }{\partial x_j}\left(u_j\left(\rho E+p\right)\right)=-\frac{\partial q_j}{\partial x_j}+\frac{\partial }{\partial x_j}\left(u_i{\tau }_{ij}\right),\hfill \end{array}$$

(1)

where $\rho$ is the fluid density, p is the pressure, E is the total energy per unit mass, $u_i$ is the ith component velocity and d is the number of geometrical dimensions. With ${\tau }_{ij}$ and $q_j$ are denoted the viscous stress tensor and the heat flux components, respectively, given by:

$$\begin{array}{c}\hfill {\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right),\\ \hfill q_j=-\frac{\mu }{Pr}\frac{\partial }{\partial x_j}\left(E+\frac{p}{\rho }-\frac{1}{2}{u}_k{u}_k\right),\end{array}$$

(2)

where $\mathsf{\mu }$ is the viscosity coefficient, $Pr$ is the Prandtl number, here equal to $0.72$, and $k=1,\dots ,d$. Finally, for a perfect gas, $p=\left(\gamma -1\right)\rho \left[E-\left(u_i{u}_i\right)/2\right],$ where $\gamma =c_p/c_v$ is the ratio of gas specific heats, here equal to $1.4$.

In compact form, the system (1) can be written as:

$$\frac{\partial \mathbf{q}}{\partial t}+\nabla ·\mathbf{F}\left(\mathbf{q},\nabla \mathbf{q}\right)=\mathbf{0},$$

(3)

where $\mathbf{q}=\{\rho ,\rho u_i,\rho E\}$ is the vector of the conservative variables, with $\mathbf{q}\in {\mathbb{R}}^m$ and $m=d+2$ and $\mathbf{F}\in {\mathbb{R}}^m\otimes {\mathbb{R}}^d$ is the sum of inviscid and viscous flux functions.

3. The Spatial Discretization: Discontinuous Galerkin Method

The discontinuous Galerkin (dG) discretization is here only briefly presented for the sake of conciseness. For a more exhaustive overview of dG methods, the interested reader is referred to [19,20].

To perform the dG discretization of Equation (3), we consider an approximation ${\mathsf{\Omega }}_h$ of the domain $\mathsf{\Omega }$, such as ${\mathsf{\Omega }}_h={\bigcup }_{K\in {\mathcal{T}}_h}K$, where ${\mathcal{T}}_h=\left\{K\right\}$ is a set of non-overlapping elements K of arbitrary shape. Moreover, we denote with ${\mathcal{F}}_h^b$ the set of faces that belong to the boundary of the discrete approximation and with ${\mathcal{F}}_h^i$ the set of internal edges of the mesh and ${\mathcal{F}}_h\stackrel{\mathrm{def}}{=}{\mathcal{F}}_h^b\cup {\mathcal{F}}_h^i$. Now, we define the solution and test function space as:

$${\mathbf{V}}_h\stackrel{\mathrm{def}}{=}\{{\mathbf{v}}_h\in L^2\left({\mathsf{\Omega }}_h\right):{\mathbf{v}}_h{|}_K\in {\mathbb{P}}^k\left(K\right)\forall K\in {\mathcal{T}}_h\},$$

(4)

where ${\mathbb{P}}^k\left(K\right)$ is the space of polynomial functions of degree at most k with no global continuity requirement, i.e., they are continuous functions only inside each element K but not on ${\mathcal{F}}_h$.

The dG formulation of Equation (3) consists in seeking ${\mathbf{q}}_h\in {\mathbf{V}}_h$ that, for an arbitrary test function ${\mathbf{v}}_h\in {\mathbf{V}}_h$, satisfies the following equation:

$$\begin{array}{cc}\hfill & {\int }_{{\mathsf{\Omega }}_h}{\mathbf{v}}_h·\frac{\partial {\mathbf{q}}_h}{\partial t}d\mathbf{x}-{\int }_{{\mathsf{\Omega }}_h}\nabla {\mathbf{v}}_h:\mathbf{F}\left({\mathbf{q}}_h,\nabla {\mathbf{q}}_h+\mathbf{r}\left(\left[\left[{\mathbf{q}}_h\right]\right]\right)\right)d\mathbf{x}\hfill \\ \hfill & +{\int }_{F\in {\mathcal{F}}_h}\left[\left[{\mathbf{v}}_h\right]\right]:\widehat{\mathbf{F}}\left({\mathbf{q}}_h^{\pm },{\left(\nabla {\mathbf{q}}_h+{\eta }_F{\mathbf{r}}_F\left(\left[\left[{\mathbf{q}}_h\right]\right]\right)\right)}^{\pm }\right)d\sigma =\mathbf{0}.\hfill \end{array}$$

(5)

In this equation, $\mathbf{r}$ and ${\mathbf{r}}_F$ are the global and local lifting operators, respectively, and ${\eta }_F$ is a stability parameter [18]. The symbol $\left[\left[·\right]\right]$ is a jump that takes into account that the functions are discontinuous on a mesh face $F\in {\mathcal{F}}_h$, and it is defined, for a generic vector $\mathbf{w}$, as:

$$\left[\left[{\mathbf{w}}_h\right]\right]\stackrel{\mathrm{def}}{=}{\mathbf{w}}_h^{+}\otimes {\mathbf{n}}_F^{+}+{\mathbf{w}}_h^{-}\otimes {\mathbf{n}}_F^{-},$$

(6)

where the ± superscripts indicate the traces of the solution over a face shared by two adjacent elements $K^{+}$ and $K^{-}$ and ${\mathbf{n}}_F^{-}$ and ${\mathbf{n}}_F^{+}$ are the unit outward normals pointing to $K^{+}$ and $K^{-}$, respectively. In the boundary integral of Equation (5), $\widehat{\mathbf{F}}$ is the sum of the inviscid and viscous numerical flux functions defined at element interfaces and boundary faces.

The inviscid numerical flux has been computed in this work using the exact Riemann solver of Gottlieb and Groth [28]. The discretization of the viscous terms is based on the BR2 scheme presented in [29] and theoretically analyzed in [18,30].

4. The Temporal Discretization

This section describes the BE-BDF2 and the RE methods and the computational aspects related to their implementation.

4.1. The BE-BDF2 Time Integration Method

Following the Method Of Lines (MOL), after that Equation (5) is numerically integrated by using the Gauss quadrature rules, it results in the following system of Ordinary Differential Equations (ODEs):

$$\mathbf{M}\frac{d\mathbf{Q}}{dt}+\mathbf{R}\left(\mathbf{Q}\right)=\mathbf{0},$$

(7)

where $\mathbf{M}$ is the global block mass matrix, $\mathbf{Q}$ is the global vector of unknown degrees of freedom and $\mathbf{R}\left(\mathbf{Q}\right)$ is the vector of residuals resulting from the integrals of the dG discretized space differential operators in Equation (5).

In this work, Equation (7) is advanced in time by using the second-order accurate BE-BDF2 scheme [11,17]. This method can be considered a variant of the one-step, two-stages TR-BDF2 scheme developed by Bank et al. [4]. In the TR-BDF2 scheme, a first solution is computed over a portion of the time-step using the Trapezoid Rule and, in the second stage, the solution at time level $t^n+\Delta t$ is computed with the BDF2 method. The BE-BDF2 method, like the TR-BDF2 scheme, is a one-step, two-stages scheme, but in the first stage the TR rule is replaced with the simpler BE method [17].

More precisely, The BE-BDF2 scheme consists in advancing Equation (7) in time from a generic time level $t^n$ to the time level $t^{n+1}=t^n+\Delta t$, successively solving the following two stages:

$$\begin{array}{c}\hfill \mathbf{M}\left({\mathbf{Q}}^{n+\gamma }-{\mathbf{Q}}^n\right)+\gamma \Delta t\mathbf{R}\left({\mathbf{Q}}^{n+\gamma }\right)=\mathbf{0},\\ \hfill \mathbf{M}\left({\mathbf{Q}}^{n+1}-\frac{1-\gamma }{\gamma }{\mathbf{Q}}^{n+\gamma }-\frac{2\gamma -1}{\gamma }{\mathbf{Q}}^n\right)+\gamma \Delta t\mathbf{R}\left({\mathbf{Q}}^{n+1}\right)=\mathbf{0}.\end{array}$$

(8)

In these equations, ${\mathbf{Q}}^n$ is a known solution at the time level $t^n$, ${\mathbf{Q}}^{n+\gamma }$ is an approximated solution computed at the time level $t^n+\gamma \Delta t$, that, for $\gamma \in \left]0,1\right[$, is a generic time level between $t^n$ and $t^{n+1}$ and ${\mathbf{Q}}^{n+1}$ is an approximated solution computed at the time level $t^{n+1}$.

The value of $\gamma$ determines the stability function of the BE-BDF2 scheme. In particular, in this work we have used $\gamma =1-\sqrt{2}/2$, a value that minimises the local truncation error of the approximation and that makes the BE-BDF2 scheme both A-stable and L-stable, i.e., unconditionally stable and stiffly accurate for stiff problems. The stability function and local truncation error of the BE-BDF2 scheme are reported in [11,17].

4.2. The Richardson Extrapolation Method

The Richardson Extrapolation (RE) is a very simple and powerful technique applicable to any numerical method [26,27,31]. More precisely, once a numerical method has been chosen for solving a system of ODEs, the application of RE results in a new numerical method with new properties, i.e., a new order of accuracy, a new local truncation error and a new stability function.

To understand the idea behind this technique it is described in its simplest form, which is the one used in this work. It can be resumed in three steps:

• Using a general time integration method of order p and a time-step size $\Delta t$, an approximation of the exact solution at time level n is computed and denoted as ${\mathbf{X}}^n$;
• Using the same time integration method but a halved time-step size, i.e., $\Delta t/2$, a better approximation of the exact solution at the same time level n is computed and denoted as ${\mathbf{Y}}^n$;
• The final approximation at time n is given by:

  $${\mathbf{Q}}^n=\frac{2^p{\mathbf{Y}}^n-{\mathbf{X}}^n}{2^p-1}.$$

  (9)

Equation (9) is just a linear interpolation of the two approximations computed at steps 1 and 2. For the BE-BDF2 scheme, since $p=2$, the RE coefficients that multiply the solutions obtained with the smaller and the larger time-step size are $4/3$ and $-1/3$, respectively.

In what follows, the “origin” of Equation (9) is briefly described, but, for a more comprehensive survey of RE, its application to several explicit and implicit methods and for the derivation of the equations reported in this subsection, please refer to [27]. Assuming that ${\mathbf{Z}}^n$ is the exact solution at time $t^n$, once that ${\mathbf{X}}^n$ and ${\mathbf{Y}}^n$ have been computed with a time integration method of order p, it is possible to write the two following equations:

$$\begin{array}{c}\hfill {\mathbf{Z}}^n-{\mathbf{X}}^n=\Delta t^p\mathbf{K}+\mathbf{O}\left(\Delta t^{p+1}\right),\\ \hfill {\mathbf{Z}}^n-{\mathbf{Y}}^n={\left(\frac{\Delta t}{2}\right)}^p\mathbf{K}+\mathbf{O}\left(\Delta t^{p+1}\right).\end{array}$$

(10)

Note that in these equations $\mathbf{K}$ depends only on the time integration method chosen and the problem to be solved, which are both the same in steps 1 and 2. The RE consists in eliminating $\mathbf{K}$ from the above two equations, thus obtaining, with very simple algebraic manipulations:

$${\mathbf{Z}}^n-\frac{2^p{\mathbf{Y}}^n-{\mathbf{X}}^n}{2^p-1}=\mathbf{O}\left(\Delta t^{p+1}\right).$$

(11)

Examining this equation, we can immediately notice that the second term on the left-hand side is exactly the approximation calculated at step 3 with Equation (9), from which it immediately follows that this approximation is of order $p+1$. We can generalize this result to iteratively use RE. In the previous example, in fact, only two numerical solutions were involved, one with the time-step $\Delta t$ and the other with the time-step $\Delta t/2$, but, if we compute other approximations with more refined grids, and we apply the same principle, an arbitrary high order accurate approximation can be found.

In the literature the RE technique is often used in two different ways: global or local extrapolation. With global extrapolation, also called passive extrapolation, the values obtained with the RE process are not used to advance the solution in time. With local extrapolation, also called active extrapolation, the values obtained with the RE process are used to calculate the next time-step. In this work, the local/active Richardson extrapolation is employed.

An important contribution to the extrapolation technique was given by Gragg [32] by considering a symmetric method, like the midpoint or the trapezoidal rules, as the “base method” for extrapolation. In this case in fact, the asymptotic error expansion of the base method contains only even powers of $\Delta t$ and the extrapolation process is particularly convenient since it eliminates two powers of $\Delta t$ at each iteration. For example, by using a symmetric method with $p=2$ and the extrapolation technique with only two grids, the final method is of order $p+2$.

The last aspect of the extrapolation technique is that it is particularly suitable for parallel implementations, since the preliminary approximations necessary to evaluate the final solution can be computed independently of each other. This last aspect is particularly attractive for large-scale problems solved by GPUs, for which an extreme parallelism in space is already employed, e.g., numerical weather prediction.

4.3. Implementation Details

In this work, since we use conservative variables and orthonormal shape functions obtained with the Gram–Schmidt process [33], the mass matrix $\mathbf{M}$, reported in Equation (8), reduces to the identity matrix $\mathbf{I}$ and therefore has been omitted in what follows. Reformulating Equation (8) in terms of a nonlinear residual function $\mathcal{F}$ and a constant vector $\mathbf{b}$, the nonlinear system of equations of a generic stage can be written in a more compact and general form as:

$$\mathcal{F}\left(\mathbf{Q}\right)=\mathbf{b},$$

(12)

with:

$$\begin{array}{cc}& \mathcal{F}\left(\mathbf{Q}\right)=\mathbf{Q}+\gamma \Delta t\mathbf{R}\left(\mathbf{Q}\right),\end{array}$$

(13)

$$\begin{array}{cc}& \mathbf{b}=-{\sum }_j{\alpha }_j{\mathbf{Q}}^j,\end{array}$$

(14)

where ${\mathbf{Q}}^j$ is a generic known solution calculated at time $t^j$ and ${\alpha }_j$ are generic stage coefficients.

Applying Newton’s scheme to Equation (12), the following linear system of equations is obtained:

$${\mathbf{A}}^i\Delta {\mathbf{Q}}^i=-\mathcal{F}\left({\mathbf{Q}}^i\right)+\mathbf{b},$$

(15)

where i is the index of the Newton iteration and ${\mathbf{A}}^i$ and $\Delta {\mathbf{Q}}^i$ are equal to:

$$\begin{array}{cc}& {\mathbf{A}}^i=\frac{\partial \mathcal{F}}{\partial \mathbf{Q}}\left({\mathbf{Q}}^i\right)=\mathbf{I}-\gamma \Delta t\frac{\partial \mathbf{R}}{\partial \mathbf{Q}}\left({\mathbf{Q}}^i\right)\end{array}$$

(16)

$$\begin{array}{cc}& \Delta {\mathbf{Q}}^i={\mathbf{Q}}^i-{\mathbf{Q}}^{i-1}.\end{array}$$

(17)

Denoting with ${\mathbf{P}}^i$ a preconditioning matrix and by applying left preconditioning to Equation (15), we solve:

$${\mathbf{P}}^i{\mathbf{A}}^i\Delta {\mathbf{Q}}^i={\mathbf{P}}^i\left(-\mathcal{F}\left({\mathbf{Q}}^i)+\mathbf{b}\right)\right).$$

(18)

In particular, at each Newton iteration, we solve Equation (18) with the restarted Generalized Minimum RESidual (GMRES) method, using a preconditioner that is based on the ILU(0) decomposition of the analytically computed ${\mathbf{A}}^i$ matrix and on the classical Block–Jacobi (BJ) approach for parallel computations. For GMRES, ILU(0) and BJ methods we rely on the algorithms available in the PETSc library [34].

To reduce the CPU time required for updating the Jacobian $\mathbf{J}=\frac{\partial \mathbf{R}}{\partial \mathbf{Q}}$ and for ILU(0) preconditioning, we compute the matrix $\mathbf{A}$ only once per time-step. This last approximation greatly contributes to speed up the solution process, since the dG solver is typically characterized by a large system matrix that is expensive to compute and factorize. However, to avoid a too slow convergence of Newton’s algorithm, the Jacobian is recomputed when $\frac{\parallel \Delta {\mathbf{Q}}^{i-1}{\parallel }_{L2}}{\parallel \Delta {\mathbf{Q}}^i{\parallel }_{L2}}\le 5$. The maximum number of Newton iterations is ${\mathrm{it}}_{\max }=10$ and the criterion of convergence of the Newton procedure is ${\epsilon }_{\mathrm{N}}={10}^{-8}$, a value chosen to minimize the temporal discretization error for any time-step size used to perform the simulations [7,35]. The GMRES parameters used are a maximum number of Krylov-subspace vectors equal to 240 with a number of restarts equal to 1 and a linear solver tolerance ${\epsilon }_{\mathrm{GMRES}}={10}^{-2}$, a value that minimizes the CPU time required by each time-step advancement [7].

Algorithms 1 and 2 provide the most important details of the BE-BDF2 and the BE-BDF2 & RE implementation, respectively. Note that, for the sake of conciseness, Algorithm 1 describes the solution process of a generic stage, since the numerical strategies adopted for both the BE-BDF2 stages are the same.

Algorithm 1 BE-BDF2

1: procedure BE-BDF2($\Delta t$, $\mathbf{A}$, ${\mathbf{Q}}^0$)  2:     $\parallel \Delta {\mathbf{Q}}^0{\parallel }_{L2}\leftarrow 1$                     ▹ Necessary value for $i=1$  3:     for $i=1,..,$ itmax do      ▹${\mathrm{it}}_{\max }$ is the maximum number of Newton iterations  4:         Solve Equation (18) to compute ${\mathbf{Q}}^i$  5:         Compute and store $\parallel \Delta {\mathbf{Q}}^i{\parallel }_{L2}$  6:         if $\parallel \Delta {\mathbf{Q}}^i{\parallel }_{L2}\le {\epsilon }_{\mathrm{N}}$ then     ▹${\epsilon }_{\mathrm{N}}$ is the Newton’s method convergence criteria  7:              Accept the solution and exit  8:          end if  9:          Compute $\eta =\frac{{\parallel \Delta {\mathbf{Q}}^{i-1}\parallel }_{L2}}{{\parallel \Delta {\mathbf{Q}}^i\parallel }_{L2}}$ 10:          if $\eta \le 5$ then 11:             $\mathbf{A}\leftarrow \mathbf{I}-\gamma \Delta t\frac{\partial \mathbf{R}}{\partial \mathbf{Q}}$      ▹ It avoids a too slow convergence of Newton’s method 12:          end if 13:     end for 14: end procedure

Algorithm 2 BE-BDF2 & RE

1: procedure BE-BDF2 and RE(${\mathbf{Q}}^n$, $\Delta t$)  ▹ Compute ${\mathbf{X}}^{n+1}$ and ${\mathbf{Y}}^{n+1}$ and then apply RE  2:     ${\mathbf{Q}}^0\leftarrow {\mathbf{Q}}^n$                    ▹ Initial guess of Newton’s method  3:     $\mathbf{A}\leftarrow \mathbf{I}-\gamma \Delta t\frac{\partial \mathbf{R}}{\partial \mathbf{Q}}$                  ▹ Compute the system matrix  4:     ${\mathbf{Q}}^{n+1}\leftarrow$ procedure BE-BDF2$\left(\Delta t,\mathbf{A},{\mathbf{Q}}^0\right)$       ▹ Compute the coarse solution  5:     ${\mathbf{X}}^{n+1}\leftarrow {\mathbf{Q}}^{n+1}$  6:     ${\mathbf{Q}}^0\leftarrow {\mathbf{Q}}^n$                    ▹ Initial guess of Newton’s method  7:     ${\mathbf{Q}}^{n+1/2}\leftarrow$ procedure BE-BDF2$\left(\Delta t/2,\mathbf{A},{\mathbf{Q}}^0\right)$    ▹ Compute an intermediate solution  8:     ${\mathbf{Q}}^0\leftarrow {\mathbf{Q}}^{n+1/2}$                ▹ Initial guess of Newton’s method  9:     ${\mathbf{Q}}^{n+1}\leftarrow$ procedure BE-BDF2$\left(\Delta t/2,\mathbf{A},{\mathbf{Q}}^0\right)$       ▹ Compute the fine solution 10:     ${\mathbf{Q}}^{n+1}\leftarrow {\mathbf{Y}}^{n+1}$ 11:     ${\mathbf{Q}}^{n+1}=\frac{4}{3}{\mathbf{Y}}^{n+1}-\frac{1}{3}{\mathbf{X}}^{n+1}$                 ▹ RE, see Equation (9) 12: end procedure

5. Numerical Results

The aim of this section is to show the behaviour and the performance of the BE-BDF2 scheme, with and without the RE technique, when employed to solve unsteady compressible inviscid and viscous flows. To this end, the following test cases have been used:

• The inviscid isentropic convecting vortex in Section 5.1;
• The laminar vortex shedding behind a circular cylinder in Section 5.2.

The isentropic vortex [9,36,37] is a well-known test case that has an exact analytical solution, here used to evaluate the accuracy and to assess the order of convergence of the investigated time integration methods. The computed order of convergence is evaluated by means of the $\eta$ error defined as:

$$\eta (\circ ,{•}_{ref})={\left({\mathsf{\Omega }}_h\right)}^{-1/2}{\parallel \circ -{•}_{ref}\parallel }_{L2},$$

(19)

where ∘ and ${•}_{ref}$ are the numerical and the "reference" solutions, respectively, with the reference solution that is set equal to the $L2$–projection of the initial solution on the dG polynomial space.

The laminar vortex shedding has been computed at a Reynolds number $Re=100$ and a Mach number $M_{\infty }=0.2$. For this test case, since an exact solution is not available, the accuracy of the time integration methods will be evaluated by comparing the results with experimental data [38] and other numerical results [39,40,41,42,43].

In order to have a fair comparison, the same Newton and GMRES parameters were used for all the implicit schemes (please see Section 4.3), and, for the explicit time integration method here used as reference, the time-step employed is the largest one allowed by the stability condition of the scheme.

5.1. The Isentropic Vortex Convection Problem

This test case consists of a uniform flow to which a perturbation is added. The initial flow field is defined by:

$$\begin{array}{ccc}\hfill u_1& =& U-\frac{\alpha M_{\infty }}{2\pi }\left(x_2-\frac{L}{2}\right)e^{\beta \left(1-r^2\right)},\hfill \\ \hfill u_2& =& U+\frac{\alpha M_{\infty }}{2\pi }\left(x_1-\frac{L}{2}\right)e^{\beta \left(1-r^2\right)},\hfill \\ \hfill T& =& 1-\frac{{\left(\alpha M_{\infty }\right)}^2\left(\gamma -1\right)}{16\beta \gamma {\pi }^2}{e}^{2\beta \left(1-r^2\right)},\hfill \\ \hfill p& =& T^{\frac{\gamma }{\gamma -1}},\hfill \end{array}$$

(20)

where T is the temperature, $U=\sqrt{\gamma }$ are the “free-stream” non dimensional velocity components and r is the distance of a generic grid point, using coordinates $\left(x_1,x_2\right)$, with respect to the vortex center, which is initially placed in the middle of the computational domain. The “free-stream” Mach number, denoted as $M_{\infty }$, is equal to $\sqrt{2}$, and $\alpha$ and $\beta$ are set equal to $5/\sqrt{2}$ and $1/2$, respectively. The computational domain is a square with side length $L=10$ that has been discretized with two uniform Cartesian grids with $25\times 25$ (coarse) and $50\times 50$ (fine) elements. Boundary conditions are periodic and the simulations are performed up to a non-dimensional final time $t_f=L$, corresponding to one period of vortex revolution using the ${\mathbb{P}}^3$ dG approximation. In what follows, the non-dimensional time-step size is defined as $\Delta t=t_f/\mathrm{ncyc}$, where ncyc is the total number of time-steps necessary to reach the end time of the simulation.

Figure 1 compares the density contours computed by several time integration schemes with the largest time-step here used, equal to $\Delta t=0.5$. In the plots the greater diffusion and dispersion errors given by BDF2 and CN2 with respect to BE-BDF2 are evident; as it is evident that, when this scheme is equipped with the RE technique, the accuracy of the solution is greatly improved.

The better accuracy of BE-BDF2 & RE is evident even when using a smaller time-step size, as highlighted in Figure 2a for $\Delta t=0.25$, while in Figure 2b, which shows the results obtained with $\Delta t=0.125$, diffusion and dispersion errors are evident only for BDF2, which is in fact the least accurate between all the schemes here reported.

A more quantitative analysis is reported in Figure 3 and

Table 1. Vortex problem—time refinement study. Simulations performed on the coarse grid using ${\mathbb{P}}^3$ dG approximation. In the table are shown density errors and the related order of convergence of different time integration methods.

Coarse Grid
$\mathbf{\Delta }\mathit{t}$	BDF2		CN2		BE-BDF2		BE-BDF2 & RE
	$\mathit{\eta }\mathbf{(}\mathit{\rho },{\mathit{\rho }}_{\mathit{ref}}\mathbf{)}$	Ord.	$\mathit{\eta }\mathbf{(}\mathit{\rho },{\mathit{\rho }}_{\mathit{ref}}\mathbf{)}$	Ord.	$\mathit{\eta }\mathbf{(}\mathit{\rho },{\mathit{\rho }}_{\mathit{ref}}\mathbf{)}$	Ord.	$\mathit{\eta }\mathbf{(}\mathit{\rho },{\mathit{\rho }}_{\mathit{ref}}\mathbf{)}$	Ord.
$5.0·{10}^{-1}$	$5.80·{10}^{-2}$	$0.47$	$5.93·{10}^{-2}$	$1.42$	$3.38·{10}^{-2}$	$1.60$	$3.01·{10}^{-3}$	$3.65$
$2.5·{10}^{-1}$	$4.20·{10}^{-2}$	$1.19$	$2.20·{10}^{-2}$	$1.85$	$1.12·{10}^{-2}$	$1.91$	$2.41·{10}^{-4}$	$3.19$
$1.25·{10}^{-1}$	$1.84·{10}^{-2}$	$1.68$	$6.11·{10}^{-3}$	$1.98$	$2.97·{10}^{-3}$	$1.99$	$2.63·{10}^{-5}$	$0.46$
$6.25·{10}^{-2}$	$5.74·{10}^{-3}$	$1.92$	$1.55·{10}^{-3}$	$2.00$	$7.48·{10}^{-4}$	$1.99$	$1.91·{10}^{-5}$	$-0.02$
$3.125·{10}^{-2}$	$1.52·{10}^{-3}$	$2.00$	$3.87·{10}^{-4}$	$1.97$	$1.88·{10}^{-4}$	$1.89$	$1.94·{10}^{-5}$	$0.00$
$1.5625·{10}^{-2}$	$3.84·{10}^{-4}$	$1.97$	$9.86·{10}^{-5}$	$1.67$	$5.07·{10}^{-5}$	$1.16$	$1.94·{10}^{-5}$	$--$
$7.8125·{10}^{-3}$	$9.81·{10}^{-5}$	$1.67$	$3.10·{10}^{-5}$	$0.61$	$2.27·{10}^{-5}$	$0.21$	$--$	$--$
$3.90625·{10}^{-3}$	$3.09·{10}^{-5}$	$0.60$	$2.03·{10}^{-5}$	$0.06$	$1.96·{10}^{-5}$	$--$	$--$	$--$
$2.00000·{10}^{-3}$	$2.04·{10}^{-5}$	$0.07$	$1.95·{10}^{-5}$	$--$	$--$	$--$	$--$	$--$
$1.00000·{10}^{-3}$	$1.95·{10}^{-5}$		$--$		$--$		$--$

and

Table 2. Vortex problem—time refinement study. Simulations performed on the fine grid using ${\mathbb{P}}^3$ dG approximation. In the table are shown density errors and the related order of convergence of different time integration methods.

Fine Grid
$\Delta \mathit{t}$	BDF2		CN2		BE-BDF2		BE-BDF2 & RE
	$\mathit{\eta }\left(\mathit{\rho },{\mathit{\rho }}_{\mathit{ref}}\right)$	Ord.	$\mathit{\eta }\left(\mathit{\rho },{\mathit{\rho }}_{\mathit{ref}}\right)$	Ord.	$\mathit{\eta }\left(\mathit{\rho },{\mathit{\rho }}_{\mathit{ref}}\right)$	Ord.	$\mathit{\eta }\left(\mathit{\rho },{\mathit{\rho }}_{\mathit{ref}}\right)$	Ord.
$5.0·{10}^{-1}$	$5.80·{10}^{-2}$	$0.47$	$5.93·{10}^{-2}$	$1.42$	$3.38·{10}^{-2}$	$1.60$	$3.03·{10}^{-3}$	$3.63$
$2.5·{10}^{-1}$	$4.20·{10}^{-2}$	$1.19$	$2.20·{10}^{-2}$	$1.85$	$1.12·{10}^{-2}$	$1.91$	$2.45·{10}^{-4}$	$3.52$
$1.25·{10}^{-1}$	$1.84·{10}^{-2}$	$1.68$	$6.11·{10}^{-3}$	$1.98$	$2.97·{10}^{-3}$	$1.99$	$2.13·{10}^{-5}$	$3.18$
$6.25·{10}^{-2}$	$5.74·{10}^{-3}$	$1.92$	$1.55·{10}^{-3}$	$2.00$	$7.48·{10}^{-4}$	$2.00$	$2.35·{10}^{-6}$	$1.42$
$3.125·{10}^{-2}$	$1.52·{10}^{-3}$	$1.98$	$3.86·{10}^{-4}$	$2.00$	$1.87·{10}^{-4}$	$2.00$	$8.78·{10}^{-7}$	$0.05$
$1.5625·{10}^{-2}$	$3.84·{10}^{-4}$	$2.00$	$9.66·{10}^{-5}$	$2.00$	$4.68·{10}^{-5}$	$2.00$	$8.49·{10}^{-7}$	$0.00$
$7.8125·{10}^{-3}$	$9.62·{10}^{-5}$	$2.00$	$2.42·{10}^{-5}$	$1.99$	$1.17·{10}^{-5}$	$1.94$	$8.50·{10}^{-7}$	$--$
$3.90625·{10}^{-3}$	$2.41·{10}^{-5}$	$1.92$	$6.10·{10}^{-6}$	$1.76$	$3.05·{10}^{-6}$	$1.46$	$--$	$--$
$2.00000·{10}^{-3}$	$6.37·{10}^{-6}$	$1.83$	$1.80·{10}^{-6}$	$--$	$1.15·{10}^{-6}$	$0.39$	$--$	$--$
$1.00000·{10}^{-3}$	$1.79·{10}^{-6}$		$--$		$8.73·{10}^{-7}$		$--$

, which summarize the results of a time refinement study performed using the coarse and the fine grid. The plots of Figure 3 report the $\eta$ errors of all conservative variables. By looking at these plots, it is evident that BE-BDF2 outperforms BDF2 and CN2. Furthermore, when this time integration scheme is equipped with the RE technique, the higher accuracy obtained with a fixed time-step, with respect to the BE-BDF2 scheme alone, is relevant, and becomes noticeable for smaller and smaller time-step sizes (for $\Delta t=1.25·{10}^{-1}$ the $\eta$ error of BE-BDF2 & RE is about two orders of magnitude smaller than the one shown by BE-BDF2). This behaviour is obviously due to the different orders of convergence of the time integration methods here investigated, whose theoretical values of 2 and 3 have been reported in the plots. Finally, note that for relatively small time-step sizes, when the spatial error overwhelms the temporal one, the $\eta$ error on each variable reaches a plateau regardless of the time scheme used, e.g., the lower bound of $\eta \left(\rho ,{\rho }_{ref}\right)$ is around $2·{10}^{-5}$ for the coarse grid and $8·{10}^{-7}$ for the fine grid.

For a more precise analysis of the order of convergence shown by the simulations, Table 1 and Table 2 report the results of the time refinement study performed with the coarse and the fine grid, respectively. In these tables have been listed the values of the $\eta$ errors and the related order of convergence only for the density variable, since the analysis for the other variables came to the same conclusions. By looking at these values we can notice that, while the schemes that are not equipped with RE reach the designed order of convergence of 2 with a monotonic increasing slope, the BE-BDF2 & RE tends towards the designed order of convergence of 3 with a monotonic decreasing slope.

5.2. The Laminar Vortex Shedding behind a Circular Cylinder

In this section are shown the results of the numerical investigation performed for the laminar flow around a two-dimensional circular cylinder with $Re=100$ and $M_{\infty }=0.2$.

The computational grid used, shown in Figure 4, is a quadratic mesh obtained with the Gmsh software [44]. It is composed of 3690 quadrilateral elements, with 72 elements lying on the cylinder surface and a wall distance of the first grid nodes around the cylinder that is approximately $5\%$ of the cylinder radius. The boundary conditions imposed are a zero-heat flux no-slip condition at the wall and the far-field boundary conditions, based on the characteristic variables, for the other boundaries. Furthermore, we have verified that the position of the far-field boundary does not significantly influence the accuracy of the numerical solution. This has been achieved by performing a numerical simulation with a computational domain that has been extended to about 700 diameters in radial direction and comparing the data obtained using this last grid with the ones shown in this section.

Starting with an already periodic initial flow field, obtained with ${\mathbb{P}}^3$ dG approximation and the explicit third order accurate SSP-RK3 scheme [45], the flow was then simulated for about 20 vortex-shedding periods using different time integration methods. In particular, the results of BE-BDF2 have been compared with respect to the ones obtained with SSP-RK3, since, due to the low Mach number of the flow, this case can be considered a mildly stiff problem for which the explicit schemes are well known to be very inefficient due to their stringent time-step restrictions. For the simulations performed with SSP-RK3, the maximum time-step allowed by stability condition has been used, i.e., about $t_{vs}/\mathrm{11,000}$ being $t_{vs}$ the vortex shedding period, while for implicit schemes the time-step size has been set as equal to $t_{vs}/50$. From these data it is evident the considerable saving of iterations that can be obtained by using the implicit schemes. Note, in fact, that the time-step size used for the BE-BDF2 and BE-BDF2 & RE schemes is about 220 times larger than the one employed for the explicit scheme.

In Figure 5 are reported the velocity magnitude contours obtained at the final simulation time. The comparison of these plots highlights that there are no significant differences between the wake flow patterns computed with the different time integration methods.

Figure 6 shows the time histories of lift and drag coefficients with a close-up of the last computed vortex-shedding period on the right plots, which demonstrate that, after about 20 vortex-shedding cycles, implicit time histories present no significant differences in phase and amplitude with respect to the corresponding explicit data, with the BE-BDF2 & RE data that are slightly closer to the explicit ones with respect to the data obtained with BE-BDF2.

A quantitative comparison is reported in

Table 3. Cylinder problem—in the table are reported the mean values of lift and drag coefficients and Strouhal numbers obtained with ${\mathbb{P}}^3$ dG approximation and different time integration methods.

	$\overline{{\mathit{C}}_{\mathit{d}}}$	$\overline{{\mathit{C}}_{\mathit{l}}}$	$\mathit{Str}$
SSP-RK3	$1.364$	$-1.18·{10}^{-4}$	$0.1615$
BE-BDF2	$1.366$	$-1.96·{10}^{-4}$	$0.1611$
BE-BDF2 & RE	$1.366$	$-3.24·{10}^{-4}$	$0.1612$

in which are listed the mean values of lift and drag coefficients and the Strouhal numbers obtained for the last 10 vortex-shedding cycles. We note that the force coefficients and the Strouhal numbers given by explicit and implicit schemes are almost the same and are in agreement with the ones obtained by other flow solvers [39,40,41,42,43]. Furthermore, the computed Strouhal values are also close to the experimental values, since, as shown in [38], the experimental measurements can be approximated with the equation $Str=0.212\left(1-\frac{21.2}{Re}\right)$. What can be concluded is that the vortex-shedding phenomenon is well-predicted using both BE-BDF2 and BE-BDF2 & RE time integration methods despite the very large time-step size here employed.

6. Conclusions

In this paper we have investigated the potential of the BE-BDF2 scheme for solving unsteady compressible flows, while also considering it in combination with the Richardson Extrapolation (RE) technique. For the isentropic vortex test case, the numerical investigation shows that BE-BDF2 outperforms in accuracy the more classical BDF2 and CN2 methods, and that, when equipped with RE, the gain in accuracy is remarkable. For the cylinder test case, considered here to evaluate the performance of BE-BDF2 and BE-BDF2 & RE for stiff-problems, the numerical investigation has shown that the vortex-shedding phenomenon and the force coefficients on the cylinder are well-predicted despite the very large time-step used, which was about 220 times larger than the one employed by the explicit SSP-RK3 scheme. These results, together with the theoretical aspects related to the BE-BDF2 formulation, have led to the following conclusions about this scheme:

• It is well suited for stiff problems thanks to its stability properties;
• Even if it is composed of two stages, the coefficient that multiplies the two Jacobian matrices is the same, i.e., $\gamma \Delta t$, therefore, in particular when it is coupled with a Matrix-Free approach and a Frozen Preconditioner strategy [8], the CPU time and memory required to advance the solution in time can be greatly reduced;
• Extrapolated values can be easily computed, exploiting the previous known solutions, and can be used as very good initial guesses for both the stages, thus accelerating the convergence of Newton’s method.

The main drawback of BE-BDF2 is that it is only second-order accurate, but this aspect can be overcome if it is used in combination with RE.

The positive aspects of BE-BDF2 & RE are:

• It is third-order accurate;
• It dramatically improves the accuracy of the solution with respect to BE-BDF2, especially for high precision requirements;
• It allows the incorporation of a strategy to adapt the time-step size, since it immediately provides error estimators evaluated on several approximation levels;
• Even if it requires the computation of two solutions, these solutions can be computed in parallel.

Encouraged by the preliminary results obtained in the present study and by the promising directions listed above, future works could contribute to the advancement of numerical simulations of unsteady flows by exploring the potential of the BE-BDF2 & RE scheme for less regular solutions and geometries, e.g., turbulent flows and real-world applications.