Abstract

This chapter presents a Finite Element solution method for the incompressible Navier- Stokes equations, in primitive variables form. To provide the necessary coupling between continuity and momentum, and enhance stability, a pressure dissipation in the form of a Laplacian is introduced into the continuity equation. The recasting of the problem variables in terms of pressure and an “auxiliary” velocity demonstrates how the effects of the pressure dissipation can be eliminated, while retaining its stabilizing properties. The method can also be interpreted as a Helmholtz decomposition of the velocity vector.

The governing equations are discretized by a Galerkin weighted residual method and, because of the modification to the continuity equation, equal interpolation for all the unknowns is permitted. Newton linearization is used and, at each iteration, the linear algebraic system is solved in a fully-coupled manner by direct or iterative solvers. For direct methods, a vector-parallel Gauss elimination method is developed that achieves execution rates exceeding 2.3 Gigaflops, i.e. over 86% of a Cray YMP-8 current peakperformance. For iterative methods, preconditioned conjugate gradient-like methods are studied and good performances, competitive with direct solvers, are achieved. Convergence of such methods being sensitive to preconditioning, a hybrid dissipation method is proposed, with the preconditioner having an artificial dissipation that is gradually lowered, but frozen at a level higher than the dissipation introduced into the physical equations.

Convergence of the Newton-Galerkin algorithm is very rapid. Results are demonstrated for two-and three-dimensional incompressible flows.

Abstract

In this article we discuss the solution of the Navier-Stokes equations modelling unsteady incompressible viscous flow, by numerical methods combining operator splitting for the time discretization and finite elements for the space discretization.

The discussion includes the description of conjugate gradient algorithms which are used to solve the advection-diffusion and Stokes type problems produced at each time step by the operator splitting methods.

Introduction and Synopsis

The main goal of this article is to review several issues associated to the numerical solution of the Navier-Stokes equations modelling incompressible viscous flow. The methodology to be discussed relies systematically on variational priciples and is definitely oriented to Galerkin approximations. Also, we shall take advantage of time discretizations by operator splitting to decouple the two main difficulties occuring in the Navier-Stokes model, namely the incompressibility condition Δ u=0 and the advection term (u.Δ)u, u being here the velocity field. The space approximation will be based on finite element methods and we shall discuss with some details the compatibility conditions existing between the velocity and pressure spaces; the practical implementation of these finite element methods will also be addressed.

This article relies heavily on [1]-[7] and does not have the pretention to cover the full field of finite element methods for the Navier-Stokes equations; concentrating on books only, pertinent references in this direction are [8]-[14] (see also the references therein).

This article is organized in sections whose list is given just below.

1. The Navier-Stokes equations for incompressible viscous flow

2. Operator splitting methods for initial value problems. Application to the Navier- Stokes equations

3. Iterative solution of the advection-diffusion sub-problems

4. […]

Summary

We describe a series of algorithms for the numerical simulation of incompressible flows. These algorithms are obtained by following a rational path from a list of design goals for practical incompressible flow solvers to their ultimate realization. Along the way, the important identity of artificial viscosity and pairs of different trial spaces for velocities and pressures is shown rigourously for the mini-element. Several numerical examples demonstrate the accuracy and versatility of the algorithms developed.

Introduction

The applications that require numerical simulations of incompressible flows may be grouped into two families:

• Engineering design and optimization: here the basic physics governing the flows to be simulated are relatively well understood, and the main requirement on the numerical methods employed is versatility, ease of use, and speed. Many configurations have to be simulated quickly, in order to develop or improve a new product. This implies that the whole process of simulating incompressible flow past an arbitrary, new configuration must take at most several days. Usually, the engineer desires a global figure, like lift and drag, as the end-product of a simulation.

• Study of basic physics’, in this case, numerical simulations are used to obtain new insight into basic physical phenomena, like vortex merging and breakdown, or the transition to turbulence. The main requirement placed on the numerical methods employed is accuracy. The geometries for which these calculations are carried out are typically very simple (boxes, channels), and the time required to perform such a simulation plays a secondary role. Some of the runs performed to date have required hundreds of CRAY-hours. Usually, the physicist desires statistical data as the end-product of such a simulation.

The Adaptive Mesh-Refinement Philosophy

In computational fluid dynamics, as well as in other problems of physics or engineering, one often encounters the difficulty that the overall accuracy of the numerical solution is deteriorated by local singularities such as, e.g., singularities near re-entrant corners, interior or boundary layers, or shocks. An obvious remedy is to refine the discretization near the critical regions, i.e., to place more grid-points where the solution is less regular. The question then is how to identify these regions automatically and how to guarantee a good balance of the number of grid-points in the refined and un-refined regions such that the overall accuracy is optimal.

Another, closely related problem is to obtain reliable estimates of the accuracy of the computed numerical solution. A priori error estimates, as provided, e.g., by the standard error analysis for finite element or finite difference methods, are in general not sufficient, since they only yield asymptotic estimates and since the constants appearing in the estimates are usually not known explicitly. Morover, they often require regularity assumptions about the solution which, for practical problems, are hardly satisfied.

Therefore, a computational fluid dynamics code should be able to give reliable estimates of the local and global error of the computed numerical solution and to monitor an automatic, self-adaptive mesh-refinement based on these error estimates.

Formation of a vortex sheet

In perfect barotropic fluid, acted upon by conservative forces with a single-valued potential, the first Helmholtz law (§1.5) says that it is not possible to endow a fluid particle with vorticity, and Kelvin's circulation theorem (§1.6) shows that the circulation around a material circuit is zero if initially zero. The question arises whether vorticity can be created without violating these theorems, and without invoking viscosity, non-conservative forces or baroclinic effects. There is no a priori reason why they are not important in subsequent motion if present initially, and so one wishes to know if vorticity can be created without appeal to these effects.

Klein [1910] addressed this question with his Kaffeelöffel experiment. (See also Betz [1950].) The conclusion is that the Helmholtz and Kelvin theorems preclude the generation of piece-wise continuous vorticity, but do not prevent the formation of vortex sheets or the generation of circulation. Consider Klein's experiment. A two-dimensional plate of width 2a is set in motion through a perfect incompressible fluid with velocity U normal to the plate. We introduce the complex potential w(z) = φ+iψ, z = x+iy. The boundary conditions are ψ = Uy on x = 0, |y| < a (the axes are taken to coincide instantaneously with the plate with the y-axis along the plate and the x-axis in the direction of motion), and w ∼ 0 as z → ∞ (circulation at infinity is not allowed).