On the accuracy of the rotation form in simulations of the Navier–Stokes equations
Introduction
The nonlinearity in the Navier–Stokes equations (NSE) can be written in several ways, which, while equivalent for the continuous NSE, lead to discretizations with different algorithmic costs, conserved quantities, and approximation accuracy, e.g. [9], [10]. These forms include the convective form, the skew-symmetric form and the rotation form, given respectively byIn turbulent flow simulations, different forms of the nonlinearity are used for different reasons. The algorithmic advantages and superior conservation properties of the rotation form (summarized in Section 2.2) have led to it being a very common choice, see, e.g. Chapter 7 in [5], [6], [24]. Herein we reveal potential limitations of rotation form, and suggest a way to overcome them.
It is known from Horiuti [14], Horiuti and Itami [15] and Zang [33] that the rotation form can lead to a less accurate approximate solution when discretized by commonly used numerical methods. Horiuti [14] and Zang [33] each give numerical experiments and accompanying analytic arguments suggesting that the accuracy loss may happen due, respectively, to discretization errors in the near wall regions (Horiuti, for finite-difference methods) and to greater aliasing errors (Zang, for spectral methods). We have noticed the same loss of accuracy in experiments from [20], [26] (for finite element methods) and suggest herein a third possibility that it is due to a combination of:
- 1.
The Bernoulli or dynamic pressure is generically much more complex than the pressure p, and thus
- 2.
Meshes upon which p is fully resolved are typically under resolved for P, and
- 3.
As the Reynolds number increases, the discrete momentum equation with either form of the nonlinearity magnifies the pressure error’s effect upon the velocity error.
We will see, for example, that in the usual formulation Velocity Error ∼ Re∗ Pressure Error, (2.13). Thus, interestingly, some of the loss of accuracy, although triggered by the nonlinear term, is due to connections between variables already present in the linear Stokes problem.
In finite element methods (FEM) the inf–sup condition for stability of the pressure places a strong condition linking velocity and pressure degrees of freedom. This condition, while quite technical when precisely stated, roughly implies that for lower order approximations the pressure degrees of freedom should correspond to the velocity degrees of freedom on a mesh one step coarser than the velocity mesh, while for higher order finite elements the polynomial degree of pressure approximations is less than the polynomial degree of velocity approximations. Thus, in either case for velocities u and Bernoulli pressures P with the same complex structures, as the mesh is refined the velocity will be often fully resolved before the Bernoulli pressure is well-resolved, see the experiments in Section 3.3 involving flow around a cylinder. Further, when an artificial problem, constructed so the kinematic pressure and Bernoulli pressure reverse complexity, is solved the observed error behavior is reversed: the convective form has much greater error than the rotation form, Section 3.2.
The question of resolution is reminiscent of Horiuti’s argument based on truncation errors in boundary layers. For example, even for a simple Prandtl-type, laminar boundary layer, the pressure p will be approximately constant in the near wall region while the Bernoulli pressure will share the boundary layer of the velocity field.
Point 3 is possibly related to aliasing errors; interestingly, the aliasing error in using different forms of the nonlinearity is governed by the resolution of the (Bernoulli or kinematic) pressure. Our suggestion of a “fix” of using grad–div stabilization (see Section 2.3) works in our tests because it addresses point 3 without requiring extra resolution.
Stabilization of grad–div type reduces the error in , see (2.15), and its (bad) scaling with respect to the Reynolds number. Moreover, since when the nonlinear terms are equivalent, this stabilization causes the three schemes to produce more closely related solutions.
Generally speaking, adding the grad–div terms to the finite element formulation is not a new idea. These terms are part of the streamline-upwinding Petrov–Galerkin method (SUPG) in [8], [12], [31]. However, in practice these terms are often omitted, and until recently it was not clear if they are needed for technical reasons of the analysis of SUPG type methods only or play an important role in computations. The role of the grad–div stabilization was again emphasized in the recent studies of the (stabilized) finite element methods for incompressible flow problems, see e.g. [2], [3], [22], [23], [25], also in conjunction with the rotation form, see [21], [26].
We shall thus compare FEM (or other variational) discretizations of the rotation form of the NSE, given bywith the usual convection form, given by:and the skew-symmetric form, given by:These are related by
Finally, we note that the rotation and convection (or skew-symmetric) forms lead to linear algebra systems with different numerical properties, which occur in time-stepping or iterative algorithms for the NSE problem. While there is an extensive literature on solvers for the convection form, see e.g. [7], not so many results are known for the rotation form. However, the few available demonstrate some interesting superior properties of the rotation form in this respect. In [26], [28] it was shown that the rotation form enables one to take into account the skew-symmetric part of the matrix in such a way that a special pressure Schur complement preconditioner is robust with respect to all problem parameters and becomes even more effective when . Such type of result is still missing for the Oseen type systems with the convective terms. An effective multigrid method for the velocity subproblem of the linearized Navier–Stokes system in the rotation form was analyzed in [27]. Finally, in [1] the special factorization of the linearized Navier–Stokes system was studied, which appears to be well suited for the rotation form.
In general, it has been reported over the years that numerical errors from a specific discretization of different forms of the nonlinear terms in the Navier–Stokes equations have different effects on the accuracy and stability of flow problems. Wilhelm and Kleiser [32] showed that for the spectral element method (in which the velocity and pressure are approximated by polynomials of order N and , respectively), numerical instabilities may occur in incompressible Navier–Stokes simulations when the rotational form is not used. Furthermore, they demonstrated that the reported instability is not caused by nonlinear effects, but it is rather a consequence of the staggered grid between velocity and pressure used in spectral element method. Analytical and numerical studies of incompressible Navier–Stokes equations of Kravchenko and Moin [19] show that aliasing errors are more destructive for spectral and high-order finite-difference calculations than for low-order finite-difference simulations of turbulent channel flow. The study shows that for skew-symmetric and rotational forms, both spectral and finite-difference methods are energy conserving even in the presence of aliasing errors. The effect on aliasing errors of the formulation of nonlinear terms for Burger’s equations and for compressible Navier–Stokes equations was examined by Blaisdell et al. [4] using a Fourier analysis and numerical experiments. The skew-symmetric form of the convective term is the form that reduced amplitude of the aliasing errors as shown theoretically and experimentally. The analysis method for the rotational form was too complicated to draw any conclusions regarding it. Alternative forms of the compressible Navier–Stokes equations were also studied from the heuristic point of view by Kennedy and Gruber [18]. Finally, we note that the name “rotation form” is sometimes attributed to the sum of two terms: [14], [15], [32], while in this paper we treat as a part of the Bernoulli pressure variable. Due to typically different discretizations (meshes) for pressure and velocity, these two alternatives lead to discrete problems with different properties. Thus, if is treated as a part of the pressure term, we expect that the under resolution of Bernoulli pressure variable may affect not only finite element, but other type of discretizations as well.
Section snippets
Differences between the nonlinearities
We now illustrate some differences between the three different forms of the NSE nonlinearity. First we discuss the Bernoulli pressure, which is used instead of usual pressure, with the rotation form of the nonlinearity, and present a bound (based on the velocity part of the Bernoulli pressure) for the rotation form FEM residual in the convective form FEM. Next, we elaborate the difference in conservation laws of the (FEM discretized) nonlinearities. Lastly in this section, we present a brief
Three numerical experiments
We consider three carefully chosen examples that, we believe, give strong support for the scenario of accuracy loss described in the introduction. We use the software FreeFem++ [13] to run the numerical tests. The models are discretized with the Crank–Nicolson method in time and with the Taylor–Hood finite elements (continuous piecewise quadratic polynomials for the velocity and linears for the pressure) in space; the nonlinear system is solved by a fixed point iteration.
Conclusions
Although the convective, skew-symmetric and rotation forms of the nonlinearity are equivalent in the continuous NSE, in finite element discretizations the rotation form may offer better physical properties (in terms of conservation laws), superior properties for iterative algorithm development, is typically more stable than the convective form, and is less expensive than computing the skew-symmetric form.
However, using rotation form requires the use of the Bernoulli pressure, which is
References (33)
- et al.
Stabilized finite element methods for the generalized Oseen problem
Comput. Methods Appl. Mech. Eng.
(2007) - et al.
The effect of the formulation of the nonlinear terms on aliasing errors in spectral methods
Appl. Numer. Math.
(1996) - et al.
Stabilized finite element methods: II. The incompressible Navier–Stokes equations
Comput. Methods Appl. Mech. Eng.
(1992) - et al.
A velocity–pressure streamline diffusion method for the incompressible Navier–Stokes equations
Comput. Methods Appl. Mech. Eng.
(1990) Comparison of conservative and rotation forms in large eddy simulation of turbulent channel flow
J. Comput. Phys.
(1987)- et al.
J. Comput. Phys.
(1998) - et al.
On the effect of numerical errors in large eddy simulations of turbulent flows
J. Comput. Phys.
(1997) A low order Galerkin finite element method for the Navier–Stokes equations of steady incompressible flow: a stabilization issue and iterative methods
Comput. Methods Appl. Mech. Eng.
(2002)- et al.
Stable and unstable formulations of the convection operator in spectral element simulations
Appl. Numer. Math.
(2000) On the rotation and skew-symmetric forms for incompressible flow simulations
Appl. Numer. Math.
(1991)
An efficient solver for the incompressible Navier–Stokes equations in rotation form
SIAM J. Sci. Comput.
Stabilized finite element schemes for incompressible flow using Scott–Vogelius elements
Appl. Numer. Math.
Spectral Methods in Fluid Dynamics
Spectral Methods Evolution to Complex Geometries and Applications to Fluid Dynamics
Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics, Numerical Mathematics and Scientific Computation
Cited by (102)
A two-level finite element method with grad-div stabilizations for the incompressible Navier–Stokes equations
2024, Journal of Computational and Applied MathematicsStabilized isogeometric collocation methods for scalar transport and incompressible fluid flow
2023, Computer Methods in Applied Mechanics and EngineeringOn augmented finite element formulations for the Navier–Stokes equations with vorticity and variable viscosity
2023, Computers and Mathematics with ApplicationsA sparse grad-div stabilized algorithm for the incompressible magnetohydrodynamics equations
2023, Computers and Mathematics with ApplicationsA parallel grad-div stabilized finite element algorithm for the Stokes equations with damping
2023, Computers and Mathematics with ApplicationsStability in 3d of a sparse grad-div approximation of the Navier-Stokes equations
2022, Journal of Mathematical Analysis and Applications
- 1
Partially supported by NSF Grant DMS 0508260.
- 2
Partially supported by the RAS program “Contemporary problems of theoretical mathematics” through the Project No. 01.2.00104588 and RFBR Grant 08-01-00415.