Introduction

Physically the Stokes equations model “slow” flows of incompressible fluids or alternatively isotropic incompressible elastic materials. In Computational Fluid Dynamics, however, the Stokes equations have become an important model problem for designing and analyzing finite element algorithms. The reason being, that some of the problems encountered when solving the full Navier-Stokes equations are already present in the more simple Stokes equations. In particular, it gives the right setting for studying the stability problem connected with the choice of finite element spaces for the velocity and the pressure. It is well known that these spaces cannot be chosen independently when the discretization is based on the “Galerkin” variational form. This method belongs to the class of saddle-point problems for which an abstract theory has been developed by Brezzi [1974] and Babuska [1973]. The theory shows that the method is optimally convergent if the finite element spaces for velocity and pressure satisfy the “Babuska-Brezzi” or “inf-sup” condition. In computations the violation of this condition often leads to unphysical pressure oscillations and a “locking” of the velocity field, cf. Hughes [1987]. During the last decade this problem has been studied thoroughly and various velocity-pressure combinations have been shown to satisfy the Babuska-Brezzi condition. Unfortunately, however, it has turned out that many seemingly natural combinations do not satisfy it. (See Girault and Raviart [1986], Brezzi and Fortin [1991], and references therein.)

In this chapter we will review a recent technique of “stabilizing” mixed methods. In this approach the standard Galerkin form is modified by the addition of mesh-dependent terms which are weighted residuals of the differential equations.

Abstract

A numerical algorithm for enforcing the conservation of mass in incompressible flow simulation is discussed and the details of the implementations in terms of standard finite volumes and finite elements are given. It is also demonstrated that both segregated and coupled iterative techniques are applicable and numerical results of test cases for two and three dimensional cavity flows are presented.

Introduction

In many applications, numerical simulations of three dimensional incompressible flows are needed. The problem of satisfying exactly the continuity equation, for these flows, is well known [1] [2].

It is conceivable to update the velocity field using the momentum equations but it is not clear how to update the pressure to conserve mass, since no pressure term appears in the continuity equation. Various methods have been introduced to tackle this problem including the penalty method, the artificial compressibility method, the artificial viscosity methods, the projection and the pressure correction methods. (The discussion, here, is limited to methods based on primitive variables. The vector potential and/or velocity vorticity formulations are not covered, see for example the work of Osswald, Ghia & Ghia [19] which still requires staggered orthogonal grids to enforce mass conservation.)

In the primitive variable methods, either staggered grids are used or the continuity equation is modified. For example in the penalty method of Temam [3], a small term proportional to the pressure is added to the continuity equation, while in Chorin's artificial compressibility method [4], the continuity equation is modified by an artificial time dependent term proportional to the time derivative of the pressure.

Numerical methods for incompressible fluid dynamics have developed to the point at which a survey of the field is both timely and appropriate. A major stimulus to the field has been the large number of applications in which incompressible flows play a crucial role, and this has spurred the interest of numerous computational engineers and mathematicians. The articles which follow provide a reasonably broad view of algorithmic and theoretical aspects of incompressible flow calculations.

It should be noted at the outset that it can be dangerous to define an algorithm for simulating incompressible flows by setting, for example, the density to be constant in a successful compressible flow algorithm. The nature of the pressure as a Lagrange multiplier rather than as a thermodynamic variable as well as the infinite speed of propagation of disturbances and other factors peculiar to incompressible flows, make algorithmic development and implementation in this context a unique undertaking (see Appendix 7A).

Perhaps the first major advance in the application of large scale digital computation to incompressible flows occured in the late 1950s with the introduction of staggered mesh techniques, exemplified, for example, by the Marker-and-Cell (MAC) scheme. The use of staggered meshes in the context of the primitive variable formulation was found to provide a stable discretization of the incompressibility constraint. Shortly thereafter, it was realized that the use of staggered meshes could be avoided by employing the streamfunction-vorticity formulation in which the incompressibility constraint does not explicitly appear. Numerous finite difference algorithms were proposed and used based on this formulation of the Navier-Stokes equations.

Introduction

The rapid development and introduction of new supercomputer systems over the last decade has opened new opportunities for numerical studies of incompressible fluid flows. A new awareness has also developed that the emerging hardware technologies influence both the nature and implementation of effective algorithms to solve these problems. Specific software implementation issues concern such varied questions as code dependencies, locality, round-off errors, storage access and capacity, input-output, workstation interfaces, cache utilization, . . . and are affected significantly by such computer hardware characteristics as the graininess of parallel systems, vector length, instruction conflicts, instruction set design, network access, distribution and paging of memory, to mention but a few. One significant result of this complex environment has been the stimulation of new ideas to make optimal use of the new supercomputer architectures and to achieve both high accuracy and high computational efficiency in the fluid simulations.

In this paper, we shall review some novel methods that are especially well suited for various aspects of incompressible fluid flow simulation studies. The key dynamical feature of these flows is the absence of shock waves, so many of the results to be stated for incompressible flows should also carry over to shock-free flows at moderate Mach numbers. Even within the context of incompressible flows, there is a wealth of dynamical phenomena to be studied, including laminar flows, transition to turbulence, turbulence, free surface flows, heat transfer, particle transport, fluid-structural interactions, and multiphase flows, among other phenomena.

Introduction

The control of fluid motions for the purpose of achieving some desired objective is crucial to many technological applications. In the past, these control problems have been addressed either through expensive experimental processes or through the introduction of significant simplifications into the analyses used in the development of control mechanisms. Only recently have flow control problems been addressed, by scientists and mathematicians, in a systematic, rigorous manner. This interest is quickly expanding so that, at this time, flow control is becoming a very active and successful area of inquiry. For example, recent publications, e.g., [1] [28], provide analyses of various aspects of flow control problems, and include one or more of the following components:

• the construction of mathematical models, invoking minimal assumptions about the physical phenomena;

• the analysis of the mathematical models to answer questions about the existence and regularity of solutions and to derive necessary conditions that optimal controls and states must satisfy;

• the construction and analysis of discretization methods for determining approximate solutions of the optimal control problems, and the rigorous derivation of error estimates; and

• the development of computer codes implementing discretization algorithms, both for the purpose of showing the efficacy of these methods, and also to solve problems of practical interest.

An optimal control or optimization problem is composed of two ingredients: a desired objective and control mechanisms that are used to (hopefully) achieve the desired objective. In a mathematical description of such problems, the desired objective is usually expressed in terms of the extremization of a functional depending on the state of the system, and possibly also on the control mechanisms.

Introduction

The finite element method has been an established method for approximating incompressible flow for a number of years. The primitive variable, velocity/pressure formulation, is the most popular way to implement the method although it is not the only possibility. Much theoretical work has been done to establish convergence and error estimates and there is a large amount of literature on the topic, see for example Temam [1984], Thomasset [1981], Girault and Raviart [1986], Gunzburger [1989]. Effort has been concentrated on two-dimensional flow and although mathematically three-dimensional flow is no more difficult, in practice the current state of computer hardware makes the implementation of three-dimensional elements much more problematic. It is only the availability of modern supercomputers that has allowed the approximation of such flows to be attempted. However, an element and method of solution that works well on a large vector processor may be quite inefficient on a fine-grained parallel computer thus the concept of the “best element” or even a “good element” may be highly dependent on the computer on which it is to be implemented.

Most methods of solution involve at least one iteration of one form or another, the innermost loop being the solution of a set of linear equations. For practical three dimensional flow problems a direct method of solution is unlikely to be a feasible proposition for almost all situations and this inner system of equations will have to be solved iteratively. “How accurately do we need to solve this system?” and “What are the interactions between the various iterations taking place?” are two questions that have to be faced by anyone implementing the finite element method for three-dimensional flow.

Introduction

This book represents a timely review of state-of-the-art algorithms in the many branches of Computational Fluid Dynamics (CFD). Algorithms for a wide range of incompressible flow phenomena are presented ranging from highly turbulent flows to low-speed non-Newtonian flows. Despite the enormous advances of the past ten years it is clear that there are still many interesting new Algorithmic directions evolving.

Much of the material presented herein is from the perspective of the algorithm researcher whose ultimate goal is the development of viable algorithms for the accurate and efficient simulation of real-world CFD problems. In this chapter, I would like to take the opportunity to discuss CFD from a different perspective - that of a commercial CFD package developer and more importantly that of the user of CFD codes. For, in the final analysis, it is the user who must apply the algorithms developed by researchers and packaged in commercial CFD codes. Over the last ten years CFD has evolved to the point where commercial CFD codes have been developed and introduced into the marketplace. These codes are being used with increasing frequency and it is probably fair to state that today the use of CFD in industry is increasingly becoming an accepted analysis tool.

Over the last decade the application of CFD has spread from its original aerospace beginnings to increasing application in a broad spectrum of industries ranging from the more traditional automotive and electronics industries to exciting new applications in the biomedical and food industries. In these industries CFD is still not used in a truly design environment; however, the real payback from the use of CFD will come when it is used on a daily basis as a design tool by design engineers who are not CFD specialists.

Introduction

Vortex methods are a type of numerical method for approximating the solution of the incompressible Euler or Navier-Stokes equations. In general, vortex methods are characterized by the following three features.

1. The underlying discretization is of the vorticity field, rather than the velocity field. Usually this discretization is Lagrangian in nature and frequently it consists of a collection of particles which carry concentrations of vorticity.

2. An approximate velocity field is recovered from the discretized vorticity field via a formula analogous to the Biot-Savart law in electromagnetism.

3. The vorticity field is then evolved in time according to this velocity field.

In the past two decades a number of different numerical methods for computing the motion of an incompressible fluid have been proposed that have the above features. In this article we consider a class of such methods which are based on the work of Chorin [1973, 1978, 1980, and 1982]. Members of this class are related by the manner in which a vorticity field in an inviscid, incompressible flow is discretized and subsequently evolved. It is common practice to use the term vortex method or the vortex method to refer to a member of this class when it is used to model the incompressible Euler equations. One can modify the vortex method and use it to model the incompressible Navier-Stokes equations by adding a random walk. This is known as the random vortex method. One can also replace the random walk by some non-random technique for solving the diffusion equation. Such methods are generally referred to as deterministic vortex methods.

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).