Computational Fluid Dynamics
The Legacy of Group T-3(taken from LA-UR-96-1426 by Norman L. Johnson)
Feel free to browse through the entire history or click on the area of interest:
The early history is presented of the prolific development of CFD methods in the Fluid Dynamics Group (T-3) at Los Alamos National Laboratory in the years from 1958 to the late 1960s. Many of the currently used numerical methods: PIC, MAC, vorticity-stream-function, ICE, ALE methods and the k-e method for turbulence' originated during this time. The rest of the paper summarizes the current research in T-3 for CFD, turbulence and solids modeling. The research areas include reactive flows, multimaterial flows, multiphase flows and flows with spatial discontinuities. Also summarized are modern particle methods and techniques developed for large scale computing on massively parallel computing platforms and distributed processors.
At times in history there often comes together a unique confluence of people and events that can change the development of history. In many ways the early development of Computational Fluid Dynamics (CFD) methods at Los Alamos National Laboratory, and in particular in the Fluid Dynamics Group, in the late 50s and through the 60s was such an example, a rare integration of unique computational resources, people and applications. Arguably, a critical factor was the creation of the worlds largest computer resources for programs of national interest that were available for exploration into alternative CFD methods. But equally, the presence of Francis H. Harlow with his prolific creativity, which continues to this day, and his colleagues were also a rare occurrence. Even though the laboratory programs at the time needed robust simulations of multimaterial, compressible flows, all applications were fair game because of the almost total absence of CFD codes at the time.
The purpose of the present work is two-fold: to review the early days of the development of CFD methods in the Fluid Dynamics group (T-3) at Los Alamos National Laboratory (at that time called Los Alamos Scientific Laboratory and hereafter referred to as Los Alamos) and to summarize the current research in T-3, as an invitation to the reader to inquire about more information as warranted.
Although this review focuses on the work from T-3, there is significant work that has been done over the years in other parts of the laboratory, some in collaboration with T-3, or, more often, as independent work. The author makes apologies for the myopic view of CFD at Los Alamos and is the first to recognize the direct contributions of colleagues in other parts of the Laboratory and the contribution of the Laboratory as a whole, in providing one of the premiere research facilities in the world.
In the early years of the Fluid Dynamics Group (T-3) in the Theoretical Division at Los Alamos, the problems of interest were multiple materials under high compression, in which solids behave like fluids. The standard approach in the 50s to numerical modeling of deforming materials was a Lagrangian treatment with staggered primary variables (thermodynamic variables at the element centers, kinematic variables at the vertices or nodes). The Lagrangian method satisfied the need for an accurate interfacial treatment, but severely suffered from mesh distortions under the large shearing deformations and instabilities. Typical simulations at the time had to be halted when the mesh entangled and painstakingly "remeshed" by hand, and then the simulation continued.
In this time of great need, the PIC method was proposed and developed by Harlow in 1957 [Evans et al.][Harlow, 1957]. The original PIC code used mass particles that carried material position, mass, and species information on a two-dimensional (2D), uniform, Eulerian mesh. It treated transient, compressible flows of multiple materials with no restrictions on interfacial deformation. It was also the first of the T-3 codes that used the technique of solution phases: the division of the computational cycle into a Lagrangian and Eulerian (remap or rezone) phase. Fig. 1-1 shows a result from an early PIC calculation. A large number of particles per element ' 16 was found to be best in 2D ' were required to reduce the inherent fluctuations of the method (this is discussed in more detail in the FLIP section 3.3 below) and, consequently, the method was memory intensive, particularly for the computers of the time (IBM 701 and 704). While the PIC scheme for fluid flow had limited application outside of Los Alamos except for plasma simulations, the T-3 PIC method did find significant use in the Soviet Union as the "Large Particle" technique. The PIC method has had a major resurgence almost three decades later in the development of FLIP (see section 3.3).
Fig. 1-1. A shock passing though two gases with a stepped interface, with a density ratio of 2 (the lighter gas is shown) . The PIC particles and mesh are not illustrated.
The success of the PIC method in solving the truly unsolvable problems of the time made the idea of forming a dedicated fluid dynamics group attractive. With the support of Stan Ulam and Conrad Longmire, the Fluid Dynamics Group (T-3) was created in the Theoretical Division in 1958, with Harlow as its first Group Leader. Harlow remained group leader until 1973. T-3 started out with seven core members, and grew to 13 members by 1963, 15 members by 1970, and 25 by 1990. Although the group was largely funded by weapons research money in these early years and weapons applications remained the main area of application, the atmosphere of the late fifties and sixties was one of free exploration of CFD techniques for solving a wide variety of applications, including incompressible, free-surface flows.
The rest of this introduction focuses on the specific techniques that were developed; these were developed with a common approach, a certain style that was characteristic to T-3. The techniques were developed under the collaboration of typically a programmer and a theorist. The involvement of a skilled programmer was essential, because each code pushed the limits of the current computer capability. The necessity for large computers precluded much of the use of this work outside of Los Alamos in the early 60s. Computer codes for each new method were written from scratch and were not intended to be exported. But, as the 60s progressed, the T-3 techniques were widely applied across the country. The development of codes for use outside Los Alamos came much later.
To address the particle fluctuations and the large memory requirements of the PIC method, the FLIC method was developed under Harlow's direction by Gentry, Martin and Daly [Gentry et al, 1966]. The FLIC method treated compressible flows of a single material on a 2D, uniform, Eulerian mesh, in which all the state variables were co-located at the cell center. The technique fluxed material across cell boundaries in the now-typical Eulerian fashion. Not surprisingly, the method suffered from stability problems from poorly coupled momentum and pressure fields, which plagued the co-located variable methods for the next three decades. The method included the capability to treat arbitrarily shaped objects by using a piece-wise linear representation of a solid boundary in the regular mesh ' a precursor to the later fractional area/volume formulation.
Fromms work was the first and only foray away from primitive variables in T-3 of velocity and pressure, and developed the first treatment of strongly contorting incompressible flows in the world: the vorticity-stream-function method for 2D, transient, incompressible flows in 1963 by Fromm and Harlow [Fromm et al, 1963]. Fromm's ideas of a "Phase-Error Correction" method anticipated monotonicity-preserving methods currently popular. The origin of this idea has been largely forgotten.
To treat incompressible, free surface flows, the MAC method was developed by Harlow and Welch [Harlow et al, 1965] as a variation of the PIC method but treating applications that extended beyond those addressed by the vorticity'stream-function method. The MAC method was the first successful technique for incompressible flows. Particles were used as markers to locate the material in the mesh and, consequently, to define the location of the free-surface. The MAC method had the advantage of a more compact finite difference stencil and tight coupling between the pressure and velocity fields. To treat the fluid incompressibility, a solution to the Poisson equation for the pressure was used. This was in contrast to later methods that solved the coupled velocity-pressure equations, as discussed by Viecelli [Viecelli, 1969]. Although the solution of Poisson's equation was numerically simple, the specification of the velocity boundary conditions were not straightforward. There was some controversy at the time about the relative stability of the MAC method, and this was resolved in the now-classic paper by Hirt [Hirt, 1968], in which he showed that the MAC method is unstable with centered momentum advection unless the viscosity is sufficiently large. This work was the precursor of the modern truncation error subtraction analysis. This controversy illustrated the T-3 approach: the development was always on the physics, with limited application of mathematical analysis of, e.g., convergence and stability properties. The MAC method is still in use and has profited from the added efficiency of modern conjugate gradient schemes for solving the Poisson equation.
Also in the late 60s, an all-speed code was developed, and called the ICE method [Amsden, 1968][Harlow et al, 1970]. The ICE method was the first approach that removed the Courant stability limitation based on the fluid sound speed. Originally, the method had a fully nonlinear implicitness, which is often replaced by a modern linear implicitness ' a more simple approach, but with the same stability properties in the limit of zero Mach number. In the limit of zero Mach number, the ICE scheme reduces to the MAC scheme.
The MAC method was the basis of the later particle-less techniques for compressible and incompressible flows embodied in the SOLA family of codes by Hirt, Nichols and others in the early 70s that were the first T-3 codes distributed internationally. These codes included extensions to two immiscible fluids in the SOLA-VOF code [Hirt et al, 1975][Hirt et al, 1981], the first broadly distributed T-3 code in its source form. One member of SOLA family, the SOLA-DF code, included a multiphase treatment with multiple velocity fields [Hirt et al, 1979][Rivard et al, 1976][Travis et al, 1979]. About the same time, the first reactive flow code, RICE, was developed by Rivard and Butler and others [Butler et al, 1976][Rivard et al, 1975], which evolved into the most widely used of the T-3 codes, APACHE-CONCHAS-KIVA lineage of codes [Ramshaw et al, 1979], discussed in the next section.
In 1967 the first 2D Lagrangian method for incompressible flows was developed in the LINC (Lagrangian-INCompressible) code. The approach was based on restricting the movement of the vertices such that the volume remained constant and, thus, was not based upon a global solution to the zero divergence of the velocity field, as in previous incompressible methods. While the LINC code was no more successful in treating flows than other multi-dimensional Lagrangian codes, its formulation led to the staggered mesh approach to coupling of the pressure and velocity fields in the MAC method and was used to explore elastic-plastic materials and surface tension effects. The most important consequence of the LINC code was the observation that mesh rezoning was needed for most problems. The second generation version of LINC, consequently, included an ALE capability. This first application of the ALE formalism paved the future for all the later ALE codes [Brackbill et al, 1973][Hirt et al, 1974], including SALE [Amsden, 1980] and its progeny. This version of LINC was also the first application of the Finite Volume method, the use of integral formulation of the conservation equations, the close cousin to the finite element methods. The Finite Volume method enabled straightforward treatment of nonorthogonal and three-dimensional meshes.
Harlow and his colleagues also contributed to the early numerical modeling of turbulence and, in particular, by the postulation of the now ubiquitous k-e model in the 60s [Harlow et al, 1967][Harlow et al, 1961]. The history of the early turbulence modeling in T-3 is included with the current modeling in section 3.5 below, in order to present a unified treatment of this complex topic.
One of the least known CFD methods outside of T-3, but one that was the precursor to the Free-Lagrangian methods, including the Smooth Particle and the Lattice Gas methods, is the PAF method, first documented in 1961 [Harlow et al, 1961]. It was the first of the mesh-less (in the sense that computational points were not associated with any mesh) and variable connectivity methods (in the sense that the connectivity changed during the simulation). It combined the lack of numerical diffusion of the Lagrangian methods with the robustness of the Eulerian methods, but without the underlying mesh and the large memory requirements of the PIC method. One way to view the method is as a molecular dynamics approach, but applied on a macroscopic scale. Computational points have a constant mass and carry all state information; they do not possess any moment of inertia, i.e., they are point masses. The particles represent parcels of fluid that interact with fluid like forces that are chosen to duplicate the equation of state and viscous effects [Harlow, 1963]. At any time, the particles interact with only their neighbors. The time evolution of the particles is just the solution of Newtons equations for a multi bodied system.
In 1965, a summary report was published [Daly et al, 1965] and comparison between fairly complex experimental data and simulations were made. The PAF technique was abandoned because of the inherent noise in the flow field as particles reconnected with different neighbors during shear flows. The PAF method also suffered from slowness of the calculation of the nearest neighbors, one which scales with N2, where N is the number of particles, if no acceleration techniques are used. Modern methods now have reduced this scaling to be linear and the approach has become computationally attractive again.
Not until 1983 were methods developed that minimized the fluctuations in the PAF method, such that even incompressible flows could be modeled [Johnson, 1983]. The smooth particle methods of recent times take a different approach and reduce this difficulty by averaging over more particles, but at the expense of less compact support and more computations. The Lattice Gas methods, which take the approach of reducing the unrestricted particle motion to moving on a regular lattice, tried to relate the fluctuations to thermal motion and averaged the solution over a large volume to eliminate the fluctuations in the macroscopic flow field.
Harlow has often said that 1968 was the last year that he could keep up with all the CFD developments around the world, so much had the entire field grown after that time. By a similar measure, the CFD methods developed in the 70s and 80s in T-3 were more application driven with close collaboration with the end users and less of explorations in CFD, an era had passed.
The full copy of this paper is "The Legacy and Future of CFD at Los Alamos" (365Kb pdf file) which is LA-UR-96-1426 by Norman L. Johnson.
See the viewgraphs for this paper. (878k, pdf format)
Get Adobe Acrobat
Mail Stop B216
|Operated by the University of California for the National Nuclear Security Administration, of the US Department of Energy. Copyright © 2001 UC | Disclaimer/Privacy|