A Strongly Consistent Finite Difference Scheme for Steady Stokes Flow and its Modified Equations

We construct and analyze a strongly consistent second-order finite difference scheme for the steady two-dimensional Stokes flow. The pressure Poisson equation is explicitly incorporated into the scheme. Our approach suggested by the first two authors is based on a combination of the finite volume method, difference elimination, and numerical integration. We make use of the techniques of the differential and difference Janet/Groebner bases. In order to prove strong consistency of the generated scheme we correlate the differential ideal generated by the polynomials in the Stokes equations with the difference ideal generated by the polynomials in the constructed difference scheme. Additionally, we compute the modified differential system of the obtained scheme and analyze the scheme's accuracy and strong consistency by considering this system. An evaluation of our scheme against the established marker-and-cell method is carried out.



There are no comments yet.



Algorithmic approach to strong consistency analysis of finite difference approximations to PDE systems

For a wide class of polynomially nonlinear systems of partial differenti...

Error analysis of the SAV-MAC scheme for the Navier-Stokes equations

An efficient numerical scheme based on the scalar auxiliary variable (SA...

Strong Consistency and Thomas Decomposition of Finite Difference Approximations to Systems of Partial Differential Equations

For a wide class of polynomially nonlinear systems of partial differenti...

Convergence and error estimates for a finite difference scheme for the multi-dimensional compressible Navier-Stokes system

We prove convergence of a finite difference approximation of the compres...

A convergent FV-FEM scheme for the stationary compressible Navier-Stokes equations

In this paper, we propose a discretization of the multi-dimensional stat...

A dynamically-consistent nonstandard finite difference scheme for the SICA model

In this work, we derive a nonstandard finite difference scheme for the S...

Numerical Simulation of Fecal Coliform Dynamics in Luruaco Lake, Colombia

The Luruaco Lake located in the Department of Atlántico, Colombia, is da...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In this paper, we consider the two-dimensional flow of an incompressible fluid described by the following system of partial differential equations (PDEs):


Here the velocities and , the pressure , and the external forces and are functions in and ; is the Reynolds number and is the Laplace operator.

A flow that is governed by these equations is denoted in the literature as a Stokes flow or a creeping flow. Correspondingly, the PDE system (1) is called a Stokes system. It approximates the Navier–Stokes system for a two-dimensional incompressible steady flow when . The last condition makes the nonlinear inertia terms in the Navier–Stokes system much smaller then the viscous forces (cf. [16], Sect. 2211), and neglecting of the nonlinear terms results in Eqs. (1). The fundamental mathematical theory of the Stokes flow is e.g. presented in [14].

Our first aim is to construct, for a uniform and orthogonal grid, a finite difference scheme for the governing system (1) which contains a discrete version of the pressure Poisson equation and whose algebraic properties are strongly consistent (or s-consistent, for brevity) [12, 9] with those of Eqs. (1). For this purpose, we use the approach proposed in [7] based on a combination of the finite volume method, numerical integration, and difference elimination. For the generated scheme we apply the algorithmic criterion to verify its s-consistency. The last criterion was designed in [12] for linear PDE systems and then generalized in [9] to polynomially nonlinear systems. The computational experiments done in papers [2, 3] with the Navier–Stokes equations demonstrated a substantial superiority in numerical behavior of s-consistent schemes over s-inconsistent ones.

The linearity of Eqs. (1) not only makes the construction and analysis of its numerical solutions much easier than in the case of the Navier–Stokes equations, but also admits a fully algorithmic generation of difference schemes for Eqs. (1) and their s-consistency verification. To perform related computations we use two Maple packages implementing the involutive algorithm (cf. [10]) for the computation of Janet and Gröbner bases: the package Janet [4] for linear differential systems and the package LDA [11] (Linear Difference Algebra) for linear difference systems.

Our second aim is to compute a modified differential system of the constructed difference scheme, i.e., modified Stokes flow, and to analyse the accuracy and consistency of the scheme via this differential system. Nowadays the method of modified equations suggested in [20] is widely used (see [6], Ch. 8 and [17], Sect. 5.5) in studying difference schemes. The method provides a natural and unified platform to study such basic properties of the scheme as order of approximation, consistency, stability, convergence, dissipativity, dispersion, and invariance. However, as far as we know, the methods for the computation of modified equations have not been extended yet to non-evolutionary PDE systems. We show how the extension can be done for our scheme by applying the technique of differential Janet/Gröbner bases.

The present paper is organized as follows. In Section 2, we generate for Eqs. (1) a difference scheme by applying the approach of paper [7]. In Section 3, we show that our scheme is s-consistent and demonstrate s-inconsistency of another scheme obtained by a tempting compactification of our scheme. The computation of a modified Stokes system for our s-consistent scheme is described in Section 4. Here, we also show by the example of the s-inconsistent scheme of Section 3 how the modified Stokes system detects the s-inconsistency. Finally, a numerical benchmark against the marker-and-cell method is presented in Section 5 and some concluding remarks are given in Section 6.

2 Difference Scheme Generation for Stokes Flow

We consider the orthogonal and uniform solution grid with the grid spacing and apply the approach of paper [7] to generate a difference scheme for Eqs. (1).

Step 1. Completion to involution (we refer to [19] and to the references therein for the theory of involution). We select the lexicographic POT (Position Over Term) [1] ranking with


Then the package Janet [4] outputs the following Janet involutive form of Eqs. (1) which is the minimal reduced differential Gröbner basis form:


We underlined the leaders, i.e., the highest ranking partial derivatives occurring in Eqs. (3). is the pressure Poisson equation which, being the integrability condition for system (1), is expressed in terms of its left-hand sides as

Remark 1

The differential polynomial in Eqs. (3) is in Eqs. (1) reduced modulo the continuity equation .

Step 2. Conversion into the integral form. We choose the following integration contour as a “control volume”

Figure 1: Integration contour (stencil ).

and rewrite equations , and into the equivalent integral form


where is the internal area of the contour .

It should be noted that we use in Eqs. (5) the original form of given in Eqs. (1) (see Remark 1) since we want to preserve at the discrete level the symmetry of system (1) under the swap transformation


Step 3. Addition of integral relations for derivatives. We add to system (5) the exact integral relations between the partial derivatives of velocities and the velocities themselves:


Step 4. Numerical evaluation of integrals. We apply the midpoint rule for the contour integration in Eqs. (5), the trapezoidal rule for the integrals (7) and approximate the double integrals as

where is the step of a square grid in the plane.

As a result, we obtain the difference equations for the grid functions

approximating functions , and the grid functions approximating partial derivatives

where :


Step 5. Difference elimination of derivatives. To eliminate the grid functions for the partial derivatives of the velocities, we construct a difference Janet/Gröbner basis form of the set of linear difference polynomials in left-hand sides of Eqs. (8) with the Maple package LDA [12] for the POT lexicographic ranking which is the difference analogue of the differential ranking used on Step 1:


The output of the LDA includes four difference polynomials not containing the grid functions . These polynomials comprise a difference scheme. Being interreduced, this scheme does not reveal a desirable discrete analogue of symmetry under the transformation (6). Because of this reason, we prefer the following redundant but symmetric form of the scheme:


where and are discrete versions of the Laplace operator acting on a grid function as

Remark 2

The difference equation of the system (10) can also be obtained (cf. [8]) from the integral form of in Eqs. (3)–(4) with the contour illustrated in Fig. 1 by using the midpoint rule for the contour integration of the and as well as for evaluation of the additional integrals


and the trapezoidal rule for the contour integration of and .

The difference polynomials (10) approximate those in Eqs. (3), and the correspondence between differential and difference Janet/Gröbner bases is a consequence of our choice of the differential (2) and difference (9) rankings.

3 Consistency Analysis

Let be the ring of differential polynomials over the field of rational functions in and . We consider the functions describing the Stokes flow (1) as differential indeterminates and their grid approximations as difference indeterminates. Respectively, we denote by the difference polynomial ring whose elements are polynomials in the grid functions with the right-shift operators and acting as translations, for example,


We denote by the differential ideal generated by the set of left-hand sides in (1) and by the difference ideal generated by the left-hand sides of Eqs. (10).

The elements in vanish on solutions of the Stokes flow (1) and those in vanish on solutions of (10). We refer to an element in (respectively, in ) as to a consequence of Eqs. (1) (respectively, of Eqs. (10)).

Definition 1

[12] We shall say that a difference equation implies the differential equation and write when the Taylor expansion about a grid point yields


It is clear that to approximate Eqs. (3), the scheme (10) must be pairwise consistent with the involutive differential form (3). We call this sort of consistency weak consistency.

Definition 2

[12] A difference polynomial set   is weakly consistent or w-consistent with differential system (3) if


The following definition establishes the consistency interrelation between the differential and difference ideals generated by Eqs. (1) and Eqs. (10), respectively. If such a consistency holds, then it provides a certain inheritance of algebraic properties of Stokes flow by the difference scheme.

Definition 3

[9] A finite difference approximation to (1) is strongly consistent or s-consistent with Stokes flow (1) if


where is a perfect difference ideal [15] generated by the elements in the difference approximation.

Theorem 3.1

[9] The s-consistency condition (17) holds if and only if a Gröbner basis of satisfies

Corollary 1

The difference scheme (10) is s-consistent with the Stokes system (1).


By its construction, the set of difference polynomials in Eqs. (10) is a Janet/Gröbner basis of the elimination ideal where is the difference ideal generated by the polynomials in Eqs. (8) (cf. [1], Thm. 2.3.4). The same set is also a Janet/Gröbner basis for the ideal and for the same POT ranking with and . It is readily verified with the LDA package. Furthermore, it is easy to see that


where are differential polynomials in Eqs. (3).

Remark 3

For the computation of the image in mapping (19) one can use the command ContinuousLimit of the package LDA.

It is clear that s-consistency of Eqs. (10) with Eqs. (1) implies w-consistency. But the converse is not true. For the numerical simulation of the Stokes flow it is tempting to replace in Eqs. (10) with a more compact discretization


Although this substitution preserves w-consistency since


the scheme is not s-consistent.

Proposition 1

The difference scheme is s-inconsistent.


The difference polynomial (20) does not belong to the difference ideal generated by the polynomial set in Eqs. (10) since is irreducible modulo the ideal . This can be shown by the direct computation of the normal form of modulo the Janet basis (10) with the routine InvReduce of the Maple package LDA.

Now let us analyse the set with respect to s-consistency. The Janet/Gröbner basis of the difference ideal computed with LDA shows that . This basis consists of seven elements. Four of them, , imply system (3) and the three remaining elements, denoted by , and , are rather cumbersome difference equations which imply, respectively, the following differential ones


and .

Equations (22) are not consequences of the Stokes equations since the differential polynomials and are irreducible modulo the differential ideal generated by the differential polynomials in Eqs. (1). It follows that there are solutions to the Stokes equations which do not satisfy Eqs. (22).

Remark 4

Equations (22) impose the limitations on the external forces which do not follow from the governing differential equations (1). This is a result of s-inconsistency.

4 Modified Stokes Flow

In the framework of the method of modified equation (cf. [17], Sect. 5.5), a numerical solution of the governing differential system (1), for given external forces and , should be considered as a set of continuous differentiable functions whose values at the grid points satisfy the difference scheme (10). Since the difference equations (10) describe the differential ones (3

) only approximately, we cannot expect that a continuous solution interpolating the grid values exactly satisfies Eqs. (

3). In reality, it satisfies another set of differential equations which we shall call the modified steady Stokes flow or modified flow for short.

Generally, the method of modified differential equation uses the representation of difference equations comprising the scheme as infinite order differential equations obtained by replacing the various shift operators in the difference equations by the Taylor series about a grid point. For equations of evolutionary type, the next step is to eliminate all derivatives with respect to the evolutionary variable of order greater than one. This step is done to obtain a kind of canonical form of the modified equation. Then, truncation of the order of the differential representations in the grid steps gives various modified equations (“differential approximations”) of the difference scheme.

As we show, the fact that both equation systems are Gröbner bases of the ideals they generate and satisfy the condition (19) of s-consistency allows to develop a constructive procedure for the computation of the modified flow. Since the finite differences in the scheme (10) approximate the partial derivatives occurring in Eqs. (3) with accuracy , it would appear reasonable that the scheme would have the second order of accuracy. For this reason, we restrict ourselves to the computation of the second order modified flow.

The Taylor expansions of the difference polynomials in Eqs. (10) at the grid point for , and at the point for read


where the terms of order are written explicitly. The calculation of the right-hand sides in Eq.(23) as well as the computation of the expressions given below was done with the use of freely available Python library SymPy (http://www.sympy.org/) for symbolic mathematics.

Remark 5

The Taylor expansions of the s-consistent difference scheme (10) and of the s-inconsistent scheme over the chosen grid points contain only the even powers of . It follows immediately from the fact that all the finite differences occurring in the equations of both schemes are the central difference approximations of the partial derivatives occurring in (3).

Furthermore, we reduce the terms of order in the right-hand sides of (23) modulo the differential Janet/Gröbner basis (10). This reduction will give us a canonical form of the second order modified flow, since given a Gröbner basis, the normal form of a polynomial modulo this basis is uniquely defined (cf. [1], Sect. 2.1). The normal form can be computed with the command InvReduce using the Maple package Janet.

Thus, the Taylor expansion of the difference polynomials yields the second order modified Stokes flow as follows:

Remark 6

Note that the symmetry under the swap transformation (6) that holds in Eqs. (23) does not hold in Eqs. (24). This symmetry breaking is a typical effect of the application of the Gröbner reduction to symmetric systems and caused by the non-symmetry of the term ordering.

As we know, Stokes flow (1) satisfies the integrability condition (4) which we rewrite as


Substitution of the Taylor expansions (24) into the equality (25) shows that the sum of the second-order terms explicitly written in formulae (24) is equal to zero. The following proposition shows that this is a consequence of the s-consistency of the scheme.

Proposition 2

Given a uniform and orthogonal solution grid with a spacing , a w-consistent difference scheme for Eqs. (3) is s-consistent only if its Taylor expansion based on the central-difference formulas for derivatives and reduced modulo system (3), after its substitution into the left-hand side of the equality (25), vanishes for every order in .


Let be a set of s-consistent difference approximations to the differential polynomials  in the Janet/Gröbner basis (3). The w-consistency of implies the central difference Taylor expansion


We consider the family of difference polynomials


with the central-difference operators , , , approximating the partial differential operators , , , with accuracy . Apparently, belongs to the perfect difference ideal generated by :

These difference operators are composed of the translations (14). For example,


with , etc., and .

From Eqs. (26) and (27), we obtain


The implication in Eq. (4) follows from the fact that the normal form of the differential polynomial (4) modulo Eqs. (3), if it is nonzero, does not belong to the differential ideal generated by the polynomials in (3) that contradicts the s-consistency of . Because of the same argument, the equality (30) holds for any .

Corollary 2

A w-consistent difference scheme for system (3) is s-consistent if and only if its set of polynomials is a difference Janet/Gröbner basis for the POT ranking (9).


” Because of our choice (9) of the ranking and the structure (3), differential Janet/Gröbner basis with the underlined leaders, a w-consistent difference scheme composed of four difference polynomials has the only difference -polynomial of the form (27) which approximates the left-hand side of the differential integrability condition (25). Together with the Taylor expansion (26), the relations (28)–(30) imply the reduction of S-polynomial (27) to zero modulo . Thus, the scheme is a Janet/Gröbner basis.

” If a w-consistent set is a Janet/Gröbner basis, then by Theorem 3.1 it is s-consistent.

We illustrate Proposition 2 and Corollary 2 by the s-inconsistent difference scheme of Section 3 where the first three difference equations coincide with those of the system (10),

and is given by Eq. (20). Because of the distinction of the last equation from in (10), the reduced Taylor expansions of equations and are different from and in system (24):