CFD-101: The Basics of Computational Fluid Dynamics Modeling
Can you imagine a computational fluid dynamics program that simulates the behavior of different materials separated by well-defined interfaces that are subject to arbitrarily large deformations? Can you also imagine this program capturing shock waves and tracking rarefactions, slip surfaces, and other non-linear hydrodynamic phenomena?
Developing such a program would be a daunting task. You may be surprised to learn that such a program was operating in 1955, long before computer graphics or mechanical pen plotters were available, and even before high-level programming languages like Fortran were popular. Fortran, or Formula Translation System, was proposed by IBM in 1954. The program having these amazing capabilities was a Particle-In-Cell (PIC) method originated by Francis H. Harlow of the Los Alamos National Laboratory (Harlow, F.H., "A Machine Calculation Method for Hydrodynamic Problems," Los Alamos Scientific Laboratory report LAMS-1956, Nov. 1955).
Central to the PIC method is the concept of a Lagrangian particle defined by a location (x,y,z). A particle is said to be Lagrangian when it moves as though it is an element of fluid. The particle may be thought of as the location of the center of mass of the fluid element. In addition to a location, Lagrangian particles are sometimes assigned one or more property values. In the PIC method, for instance, particles have specified masses and a label indicating what material they belong to.
Figure 1: PIC calculation of a 2 cm
diameter iron sphere hitting an
aluminum plate at a supersonic speed
While the underlying computational scheme used in the PIC method employs a fixed Eulerian grid, Lagrangian particles are used to move mass, momentum, and energy through this grid in a way that preserves the identities of the different materials. There are no connections between particles so they are free to move and follow the dynamics of a flow regardless of its complexity, Figure 1. Lagrangian particles are, in fact, the key feature in the PIC method that allows it to track large fluid deformations.
Why, then, isn't the PIC method more widely used for continuum fluid mechanics? For example, there are no commercial CFD programs based on this method. It could be argued that the PIC method is best for compressible flows, while most commercial applications deal with incompressible-fluid situations.
Two additional reasons why the PIC method is not more wisely used are associated with the discreteness of Lagrangian particles. It is these discrete properties and their consequences that are the subject of this note. One obvious property is that finite changes in numerical values may occur because of changes in the number of particles. The other property is less obvious and is associated with a fundamental characteristic of fluids that generally makes it difficult to track a fluid element simply by tracking its center of mass (a discrete) location.
The Discrete Problem
In the PIC method particles have finite masses. This means that when a particle moves from one control volume of the fixed Eulerian grid into another it causes discrete changes to be recorded in the mass, momentum, and energy of the cells losing and gaining the particle. Such changes introduce fluctuations in the computed values of all fluid dynamic quantities. The magnitude of the fluctuations is inversely proportional to the square root of the average number of particles in a grid cell.
Experience has shown that the PIC method works best with at least 16 particles per cell (i.e., a 4 by 4 array in two dimensions or 64 particles per cell in three dimensions). A smaller number of particles could be used when larger fluctuations could be tolerated (or when computing resources did not allow for a larger number, a frequent situation in the early days of CFD).
Experience also showed that better results were obtained when the initial placement of particles was not regular, but staggered. It is easy to see why this is so. Suppose the particles are arranged in a regular 4 x 4 array in x-y space. If the flow is only in the x direction then a column of four particles will pass from one cell to another at the same time, which would result in a very large change in the cell values. If the particles are staggered in space, however, it is more likely that only one particle at a time will cross a cell boundary, causing the minimum discrete change in cell values.
In more recent times another approach has been used to reduce the effect of discrete changes as particles move from cell to cell. This is the "smooth particle hydrodynamics" method in which particles have finite volumes that can overlap more than one grid cell at a time. As a particle approaches a cell boundary its volume continuously sweeps from one cell to the next.
Figure 2: (a) Flow in jet hitting a
wall, (b) initial particle distribution,
(c) subsequent particle distribution
showing vertical packing and
Figure 3: Same plots as in
Fig.2 B-C except the initial particle
distribution is staggered.
The Element Distortion Problem
A more difficult problem associated with Lagrangian particles is that fluid elements rarely retain simple, convex shapes. Most often a fluid element will find itself subjected to shearing, expanding, or contracting flow processes that quickly draw it out into a long ribbon-like shape. To visualize this, you might try introducing small volumes of smoke into a strong light (e.g., from a slide projector) and see how rapidly they deform into thin curtains of smoke.
This type of deformation means that material in a fluid element will not remain localized, and a Lagrangian particle following its center of mass will no longer be a good representation of the element. In a computational method element distortion can lead to a variety of problems. One of the most common problems is that particles will not retain a uniform distribution, but will tend to bunch up in some places and move apart in others.
A simple example of these processes occurs at stagnation point. Figure 2 shows what happens to a regular array of particles in a liquid jet when it strikes a wall and flows to either side of a stagnation point that is at the center of impact. The particles bunch together in a direction normal to the wall while at the same time move further apart along the wall.
If the particles in the initial distribution are staggered these deformation processes are greatly reduced. See Figure 3. Unfortunately, staggering cannot completely eliminate this problem. In other circumstances, at a separation point or in regions of strong shear, particle staggering is not sufficient to keep particles evenly distributed.
Numerical techniques can be used to add particles in expanding regions or eliminate them in regions of convergence. Or continuous repartitioning methods can be used to relocate particles for more even coverage. However, these operations introduce local smoothing that is effectively equivalent to an Eulerian computational method and throws away one of the best features of particles, namely that of their identity.
Flow separation regions cause difficulties not only because of the difficulty of maintaining a uniform particle distribution but also because of the curvature of the flow near a separation point.
To understand why flow curvature can be a problem, consider the rigid-body rotation of a fluid. Lagrangian particles placed in such a flow should move in circles about the axis of rotation. In practice this rarely happens because most particle implementations advance the location of a particle using a linear expression of velocity. For instance, the x-location of a particle at time-step n+1 would be computed as xn+1=xn+dtU, where dt is the time-step size and U is the x-component of the flow velocity at the location of the particle.
This expression, which is linear in the velocity, moves the particle in a direction tangent to the circle. Consequently, when the particle is moved along the tangent it moves to a slightly larger radius. After a sufficient number of time steps, particles will appear as though they are being thrown outward, a kind of numerical centrifugal effect.
The only way to correct for this type of behavior is to sense when the flow has curvature and to use a second-order, quadratic expression to compute new particle positions.
Diffusion processes are easy to include in particle methods using a type of random walk, or Monte Carlo model. One technique is to imagine a particle to be a point source for material that is diffusing outward. For a short time, dt, the diffusion can be represented as having a Gaussian distribution (i.e., having the solution to the diffusion equation for a point source). Since the particle cannot be subdivided, the distribution is instead treated as a probability distribution. The particle is then moved in the time interval dt to its most probable location. A random number generator is used to select a location in this probability distribution. The idea is that if enough trials are made the number of times the particle reaches a given position is proportional to the Gaussian distribution.
When particles are used as flow markers they make particularly nice graphic displays. A good example can be found in the Marker-and-Cell (MAC) method for free surface hydrodynamics (Harlow, F.H., Shannon, J.P., and Welch, J.E., "Liquid Waves by Computer," Science 149, 1092 (1965)). In this method Lagrangian particles do not carry mass but are simply used as markers to define grid regions occupied by fluid. Results produced by the MAC method have appeared in many publications to illustrate the impressive things that can be done with computational fluid dynamics.
Figure 4: Calculation of collapsing cylindrical
column of water (a) splashing over a circular
dyke (b). Particle finger looks especially
realistic, but particles were not used in computation.
Figure 4 shows a MAC-like computation of the flow of liquid originating from the collapse of a circular column (shown in outline to the left) and splashing over a cylindrical dyke. The small finger of marker particles at the top of the splash appears especially realistic. As it happens, this computation was performed using a Volume-of-Fluid (VOF) method in which Lagrangian particles had no computational role. The particles in the picture were only included in the computation to make the graphical display.
This example shows that what seems to be a strong argument for the accuracy of discrete particles, that is, their ability to capture local details, is mostly a visual effect in this case since the dynamics was computed from purely cell-averaged quantities.
Lagrangian particles are an extremely useful computational tool, especially when they are used to track small amounts of material whose dispersion is to be minimized. When particles are used as a discrete model for a continuous medium, however, it must be remembered that they have some limitations. In this sense, particles are no different than any other discrete computational method. Some of the issues that should be considered when using Lagrangian particles have been, we hope, discreetly presented in this note.