Before publication of Newton's Principia (1687), dynamics was an empirical science; i.e., it consisted of propositions that described observed behavior without any explanation for the forces that caused motion. For example, Kepler's laws are kinematic relations that describe orbital motion. Kepler (1571–1630) discovered these relations by laboriously fitting various possibilities to the voluminous measurements of planetary motion that had been recorded by Tycho Brahe. Likewise Galileo stated propositions describing the motion of freely falling bodies. The propositions are based on relating transit times for different drop heights to clear ideas of distance, time and translational velocity. These savants, however, possessed only the vaguest notion of force.

While Galileo recognized that there must be some extended cause for acceleration or retardation, he did not realize that a uniform acceleration was a consequence of a steady force; he recognized that there was a cause for acceleration but was not able to relate cause and effect. Indeed, it is difficult to imagine how there could be progress in this direction before the creation of calculus.

At the time of Galileo, the topics at the forefront of dynamics were percussion, projectile ballistics and celestial mechanics. Each of these topics had technological importance for warfare, industrial development or navigation. Percussion in particular was concerned with the terminal ballistics of musket balls as well as the effect of a forging hammer on a workpiece.

Two bodies, labeled B and B′, collide when they come together with an initial difference in velocity. Ordinarily they first touch at a point that will be termed the contact point C. During a very brief period of contact, the point C on the surface of body B is coincident with point C′ on the surface of body B′. If at least one of the bodies, B or B′, has a surface that is topologically smooth at the contact point (i.e., the surface has continuous curvature), there is a plane tangent to this surface at C; the coincident contact points C and C′ lie in this tangent plane. If both bodies are convex and the surfaces have continuous curvature near the contact point, then this tangent plane is tangential to both surfaces that touch at C; i.e., the surfaces of the colliding bodies have a common tangent plane. The direction of the normal to the tangent plane is specified by a unit vector n; this direction is termed the common normal direction. The contact force and changes in relative velocity at the contact point C will be resolved into components normal and tangential to the common tangent plane.

In practice the bodies that are colliding are composed of elastic, elastic–plastic or viscoplastic materials, so that the large contact forces acting during a collision induce both local deformations near the contact point and global deformations (vibrations) of the entire body. This chapter focuses on the local deformations in a contact region that can be represented as an elastic–perfectly plastic solid; the additional effect of global deformations will be introduced in Chapter 7.

For collisions between hard bodies, the analysis of changes in velocity during collision is simplified by assuming that the initial point of contact is surrounded by an infinitesimally small deforming region. For other purposes, however, it is necessary to consider deformations in the small region surrounding a finite area of contact.

Axial impact on a deformable body results in a disturbance which initially propagates away from the impact site at a specific speed. This disturbance is a pulse or wave of particle displacement (and consequent stress). Wave propagation relates to propagation of a coherent pulse of stress and particle displacement through a medium at a finite speed. Familiar manifestations of this phenomenon are the transmission of sound through air, water waves across the surface of the sea and seismic tremors through the earth; thus, waves exist in gases, liquids and solids. Sources of excitation may be either concentrated or distributed spatially, and brief or extended functions of time. The unifying characteristic of waves is propagation of a disturbance through a medium. Properties of the medium that result in waves and determine the speeds of propagation are the density ρ and moduli of deformability (Young's modulus E, shear modulus G, bulk modulus K, etc.).

Longitudinal Wave in Uniform Elastic Bar

Consider a uniform slender elastic bar of cross-sectional are A, elastic modulus E and density ρ; the bar contains a region with axial stress σ(x, t) that is propagating in the positive x direction as shown in Fig. 7.1.

Impact against a mechanism composed of nearly rigid bodies is a feature of many practical machines. These systems may include mechanisms where the relative velocities at joints are initially zero and finally must vanish or they can be another type of system such as a gear train or an agglomerate of unconnected bodies where at each contact the normal component of terminal relative velocity must be separating. These two classes of multibody impact problems – mechanisms and separate bodies that are touching – are distinguished from analyses of two-body impacts by the addition of constraint equations that describe limitations on relative motion at each point of contact between bodies. These constraint equations express the linkage between separate elements. Books on dynamics typically analyze the impulsive response of systems composed of rigid bodies linked by frictionless or nondissipative pinned joints.

A ball that falls in a gravitational field before colliding against a flat surface will rebound from the surface with a loss of energy that depends on the coefficient of restitution. If the ball is free, it will continue bouncing on the surface in a series of collisions; these arise because in each collision the ball is partly elastic and during the period between collisions the ball is attracted towards the surface by gravity. In Chapter 2 it was shown that an inelastic ball (0 < e* < 1) which is bouncing on a level surface in a gravitational field has both a period of time between collisions and a bounce height that asymptotically approach zero as the number of collisions increases. In other words, with increasing time this dissipative system asymptotically approaches a stable attractor – the equilibrium configuration where the ball is resting on the level surface.

Some other systems can experience energy input during each cycle of impact and flight; consequently these systems exhibit more complex behavior. For example, a pencil has a regular hexagonal cross-section with six vertices. If the pencil rolls down a plane, the mean translational speed of the axis asymptotically approaches a steady mean speed of rolling where the kinetic energy dissipated by the collision of a vertex against the plane equals the loss in gravitational potential energy as the pencil rolls from one flat side to the next.

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.