An energy-stable scheme for incompressible Navier-Stokes equations with periodically updated coefficient matrix

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

Highlights

  • Scheme incorporates the generalized Positive Auxiliary Variable (gPAV) idea and a pressure-correction type strategy.

  • Nonlinear term replaced by a linear term plus a correction term, with the latter controlled by an auxiliary variable.

  • Scheme satisfies a discrete energy stability property, irrespective of the time step size.

  • Pressure coefficient matrix is pre-computed, and velocity coefficient matrix is updated periodically.

  • Method produces accurate simulation results at large time step sizes.

Abstract

We present an energy-stable scheme for simulating the incompressible Navier-Stokes equations based on the generalized Positive Auxiliary Variable (gPAV) framework. In the gPAV-reformulated system the original nonlinear term is replaced by a linear term plus a correction term, where the correction term is put under control by an auxiliary variable. The proposed scheme incorporates a pressure-correction type strategy into the gPAV procedure, and it satisfies a discrete energy stability property. The scheme entails the computation of two copies of the velocity and pressure within a time step, by solving an individual de-coupled linear equation for each of these field variables. Upon discretization the pressure linear system involves a constant coefficient matrix that can be pre-computed, while the velocity linear system involves a coefficient matrix that is updated periodically, once every k0 time steps in the current work, where k0 is a user-specified integer. The auxiliary variable, being a scalar-valued number, is computed by a well-defined explicit formula, which guarantees the positivity of its computed values. It is observed that the current method can produce accurate simulation results at large (or fairly large) time step sizes for the incompressible Navier-Stokes equations. The impact of the periodic coefficient-matrix update on the overall cost of the method is observed to be small in typical numerical simulations. Several flow problems have been simulated to demonstrate the accuracy and performance of the method developed herein.

Introduction

This work concerns the numerical approximation of the incompressible Navier-Stokes equations in an energy-stable fashion. Energy-stable approximations are attractive in that they not only preserve the dissipative nature of the underlying continuous Navier-Stokes system, but more practically can potentially allow the use of larger time steps in computer simulations. This type of schemes are the focus of a number of previous works in the literature; see e.g. [27], [30], [32], [12], [20], [10], [25], [15], [4]. These schemes typically treat the nonlinear term fully implicitly or in a linearized fashion. Upon discretization, they would typically entail the solution of nonlinear algebraic systems within a time step, or when only a linear system needs to be solved, would involve time-dependent coefficient matrices and entail frequent re-computations (every time step) of these coefficient matrices [10]. This is a main drawback of traditional energy-stable schemes. Their computational cost per time step is typically high compared with that of semi-implicit type schemes [6], [31], [18], [16], [3], [33], [23], [14], [8], [26], which, albeit only conditionally stable, are more commonly-used in production simulations.

An interesting recent development in this area is [22], which describes a discretely energy-stable scheme employing an auxiliary energy variable in its formulation. The Navier-Stokes equations are reformulated and augmented by a dynamic equation for the auxiliary variable, which is a scalar-valued number rather than a field function. A prominent feature of this scheme lies in the reformulation of the nonlinear term,R(t)E(t)uu, where u is the velocity, R(t) is the auxiliary variable and E(t) is the shifted total kinetic energy of the system. The numerical scheme proposed in [22] treats the uu component in an explicit fashion, but controls this explicit component by an implicit treatment of R(t)E(t). The scheme is shown to satisfy a discrete energy stability property, which is also demonstrated by numerical experiments. The scheme has an interesting property that makes it computationally attractive and competitive. Within each time step it requires only the solution of linear algebraic systems with constant coefficient matrices, which can be pre-computed, for the field functions. One does need to additionally solve a nonlinear algebraic equation about a scalar-valued number. But since this nonlinear equation is about a scalar number, not a field function, its cost is very low, accounting for about a few percent of the total cost per time step [22]. A further development of this approach is discussed very recently in [21], which presents a method for treating the so-called energy-stable open boundary conditions [9], [8], [11], [7] in an energy-stable fashion on the discrete level.

While the auxiliary-variable approach and the numerical scheme from [22] possess a number of attractive properties, certain aspects of the method are less favorable and leave much to be desired. We list some of the issues here:

  • (i)

    The need for solving a nonlinear algebraic equation for the auxiliary variable is highly undesirable. While its computational cost can be negligible, the nonlinear equation causes two complications. First, the existence and uniqueness of the solution for the auxiliary variable from the discrete scheme becomes unknown. Second, the positivity of the computed values for the auxiliary variable, as is physically required by its definition, is uncertain.

  • (ii)

    The numerical scheme of [22] is formulated in a setting where the velocity and the pressure are fully coupled, and the discrete energy stability is proven in this coupled setting. When implementing the scheme, the authors have made a further approximation about the boundary vorticity, which de-couples the pressure/velocity computations in actual simulations. The stability proof, however, does not hold if this further approximation is taken into account.

  • (iii)

    It is observed in [22] that the accuracy of the method deteriorates when the time step size becomes large (or fairly large). How to improve the accuracy of the method for large (or fairly large) time step sizes, while simultaneously preserving the favorable properties that keep the computational cost relatively low, is an open issue.

  • (iv)

    In the definition of the auxiliary variable, a biased total energy (shifted by an energy constant C0) has been used in [22]. It is observed that the C0 value seems to have an influence on the accuracy of the simulation results (see the Kovasznay flow test of [22]), which is an undesirable aspect.

The first and the second issues in the above list have been addressed by [21]. In the method presented in [21], the nonlinear algebraic equation has been eliminated, and the auxiliary variable at each time step is given by an explicit formula, which ensures that its computed values are always positive. The method is formulated in a setting in which the pressure and velocity are de-coupled (barring the auxiliary variable) by a velocity-correction strategy. This scheme retains the attractive properties found in [22], such as the discrete energy stability and the need to only solve linear algebraic systems with constant pre-computable coefficient matrices.

The numerical scheme of [21] is able to achieve these important properties, in large part, thanks to the adoption of the generalized Positive Auxiliary Variable (gPAV) approach, which was originally developed in [36] for general dissipative systems. gPAV provides a means to use a general class of functions in defining the auxiliary variable, and a systematic procedure for treating dissipative partial differential equations (PDE). The gPAV procedure endows energy stability to the resultant scheme, and also can ensure the positivity of the computed values of the generalized auxiliary variable [36]. Compared with related works [22], [37], [35], [28], [34], the gPAV framework provides a more favorable way for treating the auxiliary variables, and it applies to very general dissipative systems.

In the current work we focus on the accuracy issue of the auxiliary-variable method as listed above. We would like to explore the possibility to expand its accuracy range, and aim to achieve accuracy at large (or fairly large) time step sizes, without seriously sacrificing the computational cost for incompressible Navier-Stokes equations. Summarized in this paper is our effort in this respect and a numerical scheme that largely achieves this goal.

In the current paper we present an energy-stable scheme for the incompressible Navier-Stokes equations employing the gPAV strategy. The salient feature of the scheme lies in the reformulation and numerical treatment of the nonlinear term. In the gPAV-reformulated system we replace the nonlinear term by a linear term plus a correction term, and put the correction term under control by an auxiliary variable (a scalar-valued number). Upon discretization, this leads to a velocity linear algebraic system with a coefficient matrix that can be updated periodically, in particular once every k0 time steps in the current work, where k0 is a user-specified integer parameter. The proposed scheme is observed to produce accurate results at large or fairly large time step sizes (depending on the Reynolds number). It substantially expands the accuracy range for the time step size compared with the scheme of [22] and the scheme without the current reformulation of the nonlinear term. Incidentally, we observe that this scheme is not sensitive to the energy constant C0 used in defining the auxiliary variable.

The current scheme incorporates the gPAV idea and a pressure-correction type splitting strategy, and is endowed with several attractive properties. It is energy-stable and satisfies a discrete energy stability property. No nonlinear solver is involved in this scheme, for either the field functions or the auxiliary variable. The computations for the velocity and the pressure are de-coupled. The method requires the computation of two copies of the velocity and pressure within a time step, by solving an individual de-coupled linear equation for each of them. Upon discretization, the pressure linear algebraic system involves a constant coefficient matrix that can be pre-computed, while the velocity linear system involves a coefficient matrix that can be updated periodically (every k0 time steps). The auxiliary variable is computed by a well-defined explicit formula, which guarantees the positivity of its computed values.

The coefficient-matrix update induces an extra cost, due to the re-computation and factorization involved therein. But since this is performed only occasionally (every k0 time steps), this extra cost is effectively spread over k0 time steps. In numerical simulations k0 can typically range from several dozen to several hundred, depending on the Reynolds number (k0=20 in the majority of simulations reported herein). So the impact of occasional coefficient matrix update on the overall computational cost of the current method is quite small, and can essentially be negligible in many cases.

The contribution of this work lies in the energy-stable scheme for the incompressible Navier-Stokes equations developed herein. Its favorable properties include: (i) improved accuracy, producing accurate simulation results at large (or fairly large) time step sizes; (ii) relatively low computational cost, requiring only the solution of linear systems with coefficient matrices that are pre-computable or only need to be updated periodically; and (iii) low sensitivity to the energy constant C0.

The rest of this paper is organized as follows. In Section 2 we discuss the reformulation of the incompressible Navier-Stokes equations based on the gPAV framework, and present the energy-stable scheme for the reformulated system. We prove a discrete energy stability property of the scheme, and discuss the solution algorithm and its implementation based on high-order spectral elements [29], [1], [17], [5]. In Section 3 we test the proposed method using several flow problems and investigate its accuracy, the effect of algorithmic parameters, and the computational cost. Section 4 then concludes the presentation with comments on a number of related issues.

Section snippets

Governing equations and gPAV-reformulated equivalent system

Consider some flow domain Ω (with boundary ∂Ω) in two or three dimensions, and an incompressible flow contained in Ω. The dynamics of the system is described by the incompressible Navier-Stokes equations given by, in non-dimensional form,ut+N(u)+pν2u=f,u=0, where u(x,t) is the velocity, p(x,t) is the pressure, the nonlinear term N(u)=uu, f(x,t) is an external body force, and x and t are the spatial coordinate and time. ν is the non-dimensional viscosity (inverse of Reynolds number Re),ν

Representative numerical tests

We next use several flow problems in two dimensions to test the performance of the method developed in the previous section. The spatial/temporal convergence rates of the method are first investigated using a manufactured analytic solution. Then the Kovasznay flow and the flow past a hemisphere in a narrow periodic channel are simulated to study the accuracy and stability of the method at large (or fairly large) time step sizes.

Concluding remarks

In the current paper we have developed an energy-stable scheme for simulating the incompressible Navier-Stokes equations. The scheme incorporates a pressure-correction type strategy and the generalized Positive Auxiliary Variable (gPAV) approach. The salient feature of the algorithm lies in that in the gPAV reformulated system the original nonlinear term is replaced by the sum of a linear term (M(u)) and a correction term, and the correction term is put under control by the auxiliary variable.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement

This work was partially supported by NSF (DMS-1522537).

References (37)

Cited by (7)

  • Stable and decoupled schemes for an electrohydrodynamics model

    2023, Mathematics and Computers in Simulation
    Citation Excerpt :

    Recall that the auxiliary variable approach has been found to be a fairly general tool to propose efficient schemes for some types of equations. It has first successfully used for gradient flows [5,7,16,18,19,27,29], then extended to the NS equations or other dissipative or conservative systems; see, e.g., [8,9,9–11,28]. Our main contribution is the construction and analysis of a class of decoupled stable schemes is based on the reformulation of the electrohydrodynamic system using auxiliary functions.

  • An unconditionally energy-stable scheme for the convective heat transfer equation

    2023, International Journal of Numerical Methods for Heat and Fluid Flow
View all citing articles on Scopus
View full text