The simulation of phenomena governed by the two-dimensional Euler equations are the first and simplest example in which vortex methods have been successfully used. The reason can be found in Kelvin's theorem, which states that the circulation in material – or Lagrangian – elements is conserved. Mathematically, this comes from the conservative form of the vorticity equation. Following markers – or particles – where the local circulation is concentrated is thus rather natural. At the same time, the nonlinear coupling in the equations resulting from the velocity evaluation immediately poses the problem of the mollification of the particles into blobs and of the overlapping of the blobs, which soon was realized to be a central issue in vortex methods.

The two-dimensional case thus encompasses some of the most important features of vortex methods. We first introduce in Section 2.1 the properties of vortex methods by considering the classical problem of the evolution of a vortex sheet. We present in particular the results obtained by Krasny in 1986 [129, 130] that demonstrated the capabilities of vortex methods and played an important role in the modern developments of the method. We then give in Sections 2.2 to 2.4 a more conventional exposition of vortex methods and of the ingredients needed for their implementation: choice of cutoff functions, initialization procedures, and treatment of periodic boundary conditions. Section 2.6 is devoted to the convergence analysis of the method and to a review of its conservation properties.

The goal of this book is to present and analyze vortex methods as a tool for the direct numerical simulation of incompressible viscous flows. Its intended audience is scientists working in the areas of numerical analysis and fluid mechanics. Our hope is that this book may serve both communities as a reference monograph and as a textbook in a course of computational fluid dynamics in the schools of applied mathematics and engineering.

Vortex methods are based on the discretization of the vorticity field and the Lagrangian description of the governing equations that, when solved, determine the evolution of the computational elements. Classical vortex methods enjoy advantages such as the use of computational elements only in cases in which the vorticity field is nonzero, the automatic adaptivity of the computational elements, and the rigorous treatment of boundary conditions at infinity. Until recently, disadvantages such as the computational cost and the inability to treat accurately viscous effects had limited their application to modeling the evolution of the vorticity field of unsteady high Reynolds number flows with a few tens to a few thousands computational elements. These difficulties have been overcome with the advent of fast summation algorithms that have optimized the computational cost and recent developments in numerical analysis that allow for the accurate treatment of viscous effects. Vortex methods have reached today a level of maturity, offering an interesting alternative to finite-difference and spectral methods for high-resolution numerical solutions of the Navier–Stokes equations.

The need of a specific discussion of vortex schemes in the context of three-dimensional flows stems from the very nature of the vorticity equation that in three dimensions incorporate a stretching term. This term fundamentally affects the dynamics of the flow; it is in particular responsible for vorticity intensification mechanisms that make long-time inviscid calculations very difficult. Vorticity stretching is considered as the mechanism by which energy is being transferred between the large and the small scales in the flow. In order to resolve related phenomena, such as the energy cascade, an adequate treatment of diffusion is thus even more crucial than in two dimensions. However, the recipes for deriving diffusion algorithms are the same in two and three dimensions (they are discussed in Chapter 5), and we focus here on inviscid three-dimensional vortex schemes. Vorticity intensification in general is associated with a rapid stretching of Lagrangian elements, which makes it also crucial to maintain the regularity of the particle mesh; we refer to Chapter 7 for a general discussion of regridding techniques.

We will discuss here two classes of vortex methods that extend to three dimensions the two-dimensional schemes introduced in Chapter 2. In the first one, the vorticity is replaced by a set of points (particles), just as in two dimensions, but these particles carry vectors instead of scalars. The stretching term in the vorticity equation is accounted for by appropriate laws that modify the circulations of the particles. We call these methods vortex particle methods.

In this chapter we present boundary conditions for the vorticity–velocity formulation of the Navier–Stokes equations and we describe their implementation in the context of vortex methods. We restrict our discussion to flows bounded by impermeable, solid walls, although several of the ideas can be extended to other cases such as free-surface flows.

The direct numerical simulation of wall-bounded flows requires accurately resolving the unsteady physical processes of vorticity creation and evolution in small regions near the boundary. Vortex methods directly resolve the vorticity field, and they automatically adapt to resolve strong vorticity gradients in regions near the wall, but they are faced with the algorithmic complication of dealing with the no-slip boundary condition. The no-slip boundary condition is expressed in terms of the velocity field at the wall and does not involve explicitly the vorticity.

Mathematically we may understand this difficulty by considering the kinematic and dynamic description of the flow motion and observing that there is an inconsistency between the number of equations and the number of boundary conditions. The kinematic description of the flow, relating the velocity to the vorticity, is an overdetermined set of equations if we prescribe all the components of the velocity at the boundary. On the other hand, no vorticity boundary condition is readily available for the Navier–Stokes equations that govern the dynamic description of the flow.

Physically, the no-slip boundary condition expresses the requirement that the flow field must adhere to the boundary.

Vortex methods were initially conceived as a tool to simulate the inviscid dynamics of vortical flows. The vorticity carried by the fluid elements is conserved in inviscid flows and simulating the flow amounts to the computation of the velocity field. In bounded domains the velocity field is constrained by the conditions imposed by the type and the motions of the boundaries. For an inviscid flow, it is not possible to enforce boundary conditions for all three velocity components as we have lost the highest-order viscous term from the set of governing Navier–Stokes equations. Usually for inviscid flows past solid bodies we impose conditions on the velocity component locally normal to the boundary.

The description of an inviscid flow can be facilitated when the velocity field is decomposed into two components that have a kinematic significance. In this decomposition, a rotational component accounts for the velocity field due to the vorticity in the flow whereas a potential component is used in order to enforce the boundary conditions and to ensure the compatibility of the velocity and the vorticity field in the presence of boundaries. This is the well-known Helmholtz decomposition.

Alternatively the evolution of the inviscid flow can be described in terms of an extended vorticity field. The enforcement of a boundary condition for the velocity components normal to the boundary does not constrain the wall-parallel velocity components. This allows for velocity discontinuities across the interface that may be viewed as velocity gradients over an infinitesimal region across the boundary.

In vortex methods the flow field is recovered at every location of the domain when one considers the collective behavior of all computational elements. The length scales of the flow quantities that are been resolved are characterized by the particle core rather than the interparticle distance. These observations, which stem from the definition itself of vortex methods and are confirmed by its numerical analysis, differentiate particle methods from schemes such as finite differences.

The essense of the method is the “communication” of information between the particles, that requires a particle overlap. As a result, a computation is bound to become inaccurate once the particles cease to overlap. Computations involving nonoverlapping finite core particles should be regarded then as modeling and not as direct numerical simulations. Excluding case-specific initial particle distributions (e.g., particles placed on concentric rings to represent an azimuthally invariant vorticity distribution) the loss of overlap (and excessive overlap) is an inherent problem of purely Lagrangian methods.

The cause of the problem is the flow strain that may cluster particles in one direction and spread them in another in the neighborhood of hyperbolic points of the flow map, resulting in nonuniform distributions. At the onset of such particle distributions no error is usually manifested in the global quantities of the flow such as the linear and the angular impulse. However, locally the vorticity field becomes distorted and spreading of the particles results in loss of naturally present vortical structures, whereas particle clustering results in the appearance of unphysical ones on the scale of the interparticle separation.

A fundamental issue in the use of vortex methods is the ability to use efficiently large numbers of computational elements for simulations of viscous and inviscid flows.

The traditional cost of the method scales as O(N2) as the N computational elements and particles induce velocities at each other, making the method unacceptable for simulations involving more than a few tens of thousands of particles. We reduce the computation cost of the method by making the observation that the effect of a cluster of particles at a certain distance may be approximated by a finite series expansion. When the space is subdivided in uniform boxes it is straightforward to construct an O(N3/2) algorithm [189]. In the past decade faster methods have been developed that have operation counts of O(N log N) [17] or O(N) [91], depending on the details of the algorithm. In these algorithms the particle population is decomposed spatially into clusters of particles (see, for example, Figure B.1) and we build a hierarchy of clusters (a tree data structure) – smaller neighboring clusters combine to form a cluster of the next size up in the hierarchy and so on. The hierarchy allows one to determine efficiently where the multipole approximation of a certain cluster is valid.

The N-body problem appears in many fields of engineering and science ranging from astrophysics to micromagnetics and computer animation. In the past few years these N-body solvers have been implemented and applied in simulations involving vortex methods.

The goal of this appendix is to provide the mathematical background needed in the numerical analysis of vortex methods and, more generally, particle methods.

The two first sections are devoted to particle approximations of the solutions to advection equations. The third section summarizes some mathematical features of the Navier–Stokes equations. The numerical analysis carried out in Chapters 2 and 3 results from a combination of the results hereafter derived.

As we have seen in several occasions, the basic feature of vortex methods is that the data, that is the initial vorticity and the source terms at the boundary, are discretized on Lagrangian elements where the circulation is concentrated. These elements are termed particles, and mathematically they consist in delta functions. In Section A.1 we answer the following question: in which sense can a set of particles be used to approximate a given smooth function? We then proceed in Section A.2 to demonstrate that particles moving along a given flow are explicit exact weak solutions, in a sense that we precisely define, of the corresponding advection equation. This is the mathematical reason why particle methods are suitable for the numerical approximation of transport equations. We then show stability properties for the weak form of the transport equation. Together with the results of Section A.1, these stability estimates are the central tool for the numerical analysis of particle methods for linear equations.

In numerical simulations it is desirable to use numerical methods that are well suited to the physics of the problem at hand. As the dominant physics of a flow can vary in different parts of the domain, it is often advantageous to implement hybrid numerical schemes.

In this chapter we discuss hybrid numerical methods that combine, to various extents, vortex methods with Eulerian grid-based schemes. In these hybrid schemes, Lagrangian vortex methods and Eulerian schemes may be combined in the same part of the domain, in which each method is used in order to discretize different parts of the governing equations. Alternatively, vortex methods and grid-based methods can be combined in the same flow solver, in which each scheme resolves different parts of the domain. In this case we will discuss domain-decomposition formulations. Finally we consider the case of using different formulations of the governing equations in different parts of the domain. In that context we discuss the combination of the velocity–pressure formulation (along with grid-based methods) and the velocity–vorticity formulation (along with vortex methods) for the governing Navier–Stokes equations.

For simplicity, we often use in this chapter the terminology of finite-difference methods but it must be clear that in most cases the ideas can readily be extended to other Eulerian methods, such as finite-element or spectral methods.

One of the attractive features of vortex methods is the replacement of the nonlinear advection terms with a set of ordinary differential equations for the trajectories of the Lagrangian elements, resulting in robust schemes with minimal numerical dissipation.

Vortex methods were originally conceived as a tool to model the evolution of unsteady, incompressible, high Reynolds number flows of engineering interest. Examples include bluff-body flows and turbulent mixing layers. Vortex methods simulate flows of this type by discretizing only the vorticity-carrying regions and tracking the computational elements in a Lagrangian frame. They provide automatic grid adaptivity and devote little computational effort to regions devoid of vorticity. Moreover the particle treatment of the convective terms is free of numerical dissipation.

Thirty years ago simulations using inviscid vortex methods predicted the linear growth in the mixing layer and were able to predict the Strouhal frequency in a variety of bluff-body flow simulations. In three dimensions, we have seen that inviscid calculations using the method of vortex filaments have provided us with insight into the evolution of jet and wake flows. However, the inviscid approximation of high Reynolds number flows has its limitations. In bluff-body flows viscous effects are responsible for the generation of vorticity at the boundaries, and a consistent approximation of viscous effects, including diffusion, is necessary at least in the neighborhood of the body. In three-dimensional flows, vortex stretching and the resultant transfer of energy to small scales produce complex patterns of vortex lines. The complexity increases with time, and viscous effects provide the only limit in the increase of complexity and the appropriate mechanism for energy dissipation. In this chapter we discuss the simulation of diffusion effects in the context of vortex methods.