The earliest successful application of computational fluid dynamics (CFD) was in connection with the Manhattan Project during World War II. Researchers used computations to study the propagation and interaction of shock waves, a subject crucial to the success of the atomic bomb.

Shock Wave Discontinuities

Shock waves are mathematically treated as discontinuities, but it was quickly recognized this would cause problems for any numerical solution. Of course, a shock wave is not a true physical discontinuity, but a very narrow transition zone whose thickness is on the order of a few molecular mean-free paths. Application of the conservation of mass, momentum, and energy conditions across a shock wave requires that there be a transformation of kinetic energy into heat energy. Physically, this transformation can be represented as a viscous dissipation, which was the idea latched onto by early investigators.

By introducing an unphysically large value of viscosity, investigators were able to thicken shock transition zones to where they could be resolved computationally. This artificial increase in the value of viscosity became known as an artificial viscosity.

Viscosity Dissipation

If the viscosity isn’t large enough, velocity oscillations about the correct mean velocity are observed to develop behind a shock. These oscillations can be interpreted as a macroscopic version of heat energy, i.e., fluctuating kinetic energy in place of fluctuating molecular energy. In hydraulic jumps, the hydraulic analogy of a shock wave, this fluctuating energy appears as a sequence of large eddies behind the jump.

Artificial Viscosity

The proper formulation and magnitude needed for an artificial viscosity has undergone many refinements over the years and includes tests to apply this viscosity only in regions undergoing strong compression and with magnitudes that are various functions of the first and/or second power of the compression rate. The culmination of these refinements is best exhibited in the method pioneered by Godunov (S.K. Godunov, Mat. Sbornik 47, 271 (1959); translated as JPRS 7225, U.S. Dept. Com., 1960), in which a local “shock tube” or elementary wave solution is used to capture the existence and propagation characteristics of shock and rarefaction waves.

Numerical Viscosity

Although artificial viscosity was introduced for numerical reasons, it is an elective addition used to modify a physical process so that it can be more easily computed. Artificial viscosity should not be confused with numerical viscosity, which is an unwanted consequence of certain types of numerical approximations.

Numerical viscosity arises from discrete approximations to the momentum advection terms in Eulerian equations, or from re-zoning operations used in Lagrangian formulations. The origin of the effect is the use of a homogenizing assumption in the elements or control volumes underlying the approximation scheme. For instance, when momentum is exchanged between neighboring elements through a convective flux the resulting contributions in a given element are combined with the momentum already there to arrive at a new value of average momentum for that element. This combining or homogenizing process introduces a smoothing effect. When another step to advance time is taken, this new value is passed on to the next element in the direction of flow. Repetition of this smoothing operation over the many steps needed to carry a solution forward in time contributes to a “diffusion” of momentum in the direction of flow.

Strictly speaking, the numerical diffusion does not behave like a true viscous diffusion because it is associated with fluid convection and does not possess the correct stress-versus-strain-rate dependency associated with a real viscosity. For example, numerical diffusion does not satisfy Newtonian relativity because it depends on the choice of computational grid, which is an absolute reference frame for the numerical approximations. Also, because the amount of numerical diffusion is proportional to the velocity of flow through a grid, it does not have the rotational symmetry possessed by a real viscosity.

Numerical Approximations

Research into numerical approximation schemes that minimize numerical viscosity effects is a continuing activity of a large part of the CFD community. The difficulty in developing such schemes is that some smoothing must always be incorporated into a numerical solution to keep it computationally stable and to smooth out dispersion errors. Dispersion errors are those errors that arise because components of a solution having different grid resolution requirements may propagate through the grid with slightly different speeds. Whenever this occurs, unphysical oscillations develop in the solution where these components reinforce or cancel one another.

The trick is to develop approximation schemes that remain accurate (i.e., have a minimum of numerical smoothing) and at the same time are robust (i.e., have sufficient numerical smoothing to compensate for dispersion errors and to remain computationally stable for a wide range of applications).

What FLOW-3D Does

In FLOW-3D the default method is a first-order, upstream, advection technique that is extremely robust, but which introduces some numerical viscosity. If it is determined that this numerical viscosity is excessive, because sharp velocity profiles must be computed without the luxury of high grid resolution, then a second-order accurate, monotonicity preserving option can be employed with the flip of a switch.

For compressible flows, an implicitly coupled pressure-velocity solution option can be used in FLOW-3D to capture shock waves and minimize the appearance of post-shock oscillations.