The Legacy and Future of CFD At Los Alamos

Summary

This paper reviews the prolific early history and modern developments of Computational Fluid Dynamics (CFD) methods created in the Fluid Dynamics Group (T-3) at Los Alamos National Laboratory. It details foundational methods developed from 1958 through the late 1960s—including PIC, MAC, ICE, ALE, and k-ε turbulence models—and surveys modern research areas such as reactive combustion modeling (KIVA), multiphase flows (CFDLIB), particle methods (FLIP), casting simulation (TELLURIDE), and parallel computing for global climate and mantle dynamics.

Cover Page (LA-UR-96-1426)

LA-UR-96-1426- Approved for public release; distribution is unlimited.

Title: The Legacy and Future of CFD At Los Alamos Author(s): Norman L. Johnson Published in: Proceedings of the 1996 Canadian CFD Conference Invited talk. Ottawa, Canada June 3-4, 1996

Los Alamos NATIONAL LABORATORY Los Alamos National Laboratory, an affirmative action/equal opportunity employer, is operated by the University of California for the U.S. Department of Energy under contract W-7405-ENG-36. By acceptance of this article, the publisher recognizes that the U.S. Government retains a nonexclusive, royalty-free license to publish or reproduce the published form of this contribution, or to allow others to do so, for U.S. Government purposes. Los Alamos National Laboratory requests that the publisher identify this article as work performed under the auspices of the U.S. Department of Energy. Los Alamos National Laboratory strongly supports academic freedom and a researcher’s right to publish; as an institution, however, the Laboratory does not endorse the viewpoint of a publication or guarantee its technical correctness. Form 836 (10/96)

Page 2

Intentionally blank

Abstract and Introduction

THE LEGACY AND FUTURE OF CFD AT LOS ALAMOS Norman L. Johnson Group T-3, Theoretical Division Los Alamos National Laboratory Los Alamos, NM 87545 [email protected]

Abstract 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 1960’s. Many of the currently used numerical methods –PIC, MAC, vorticity-stream-function, ICE, ALE methods and the k-ε 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.

1. Introduction 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 50’s and through the 60’s was such an example – a rare integration of unique computational resources, people and applications. Arguably, a critical factor was the creation of the world’s 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. The first resource for further information is to explore our Web page (http://gnarly.lanl.gov/home.html), which provides contact information for each area of application.

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.

2. Early History of Group T-3 (1958-68) 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 50’s 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.

PIC: Particle-In-Cell In this time of great need, the PIC method was proposed and developed by Harlow in 1957 [42][46]. 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

Early History: PIC, FLIC, Vorticity-Stream Function, and MAC

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) [47]. The PIC particles and mesh are not illustrated.]

Formation and Style of Group T-3 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 60’s. Computer codes for each new method were written from scratch and were not intended to be exported. But, as the 60’s progressed, the T-3 techniques were widely applied across the country. The development of codes for use outside Los Alamos came much later.

FLIC: Fluid-In-Cell 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 [45]. 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.

Vorticity and Stream Function Method Fromm’s 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 [44]. Fromm’s ideas of a “Phase-Error Correction” method anticipated monotonicity-preserving methods currently popular. The origin of this idea has been largely forgotten.

MAC: Marker-And-Cell To treat incompressible, free surface flows, the MAC method was developed by Harlow and Welch [50] as a variation of the PIC method but treating applications that extended beyond those addressed by the vorticity–

Early History: ICE, SOLA, LINC, ALE, Turbulence, PAF

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 [97]. 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 [58], 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.

ICE: Implicit-Continuous-Fluid-Eulerian Also in the late 60’s, an “all-speed” code was developed, and called the ICE method [8][53]. 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.

SOLA and Reactive Flow Codes 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 70’s that were the first T-3 codes distributed internationally. These codes included extensions to two immiscible fluids in the SOLA-VOF code [55][57], 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 [56][86][96]. About the same time, the first reactive flow code, RICE, was developed by Rivard and Butler and others [27][87], which evolved into the most widely used of the T-3 codes, APACHE-CONCHAS-KIVA lineage of codes [85], discussed in the next section.

LINC (Lagrangian-Incompressible) and ALE (Arbitrary-Lagrangian-Eulerian) 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 [16] [54], including SALE [9] 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.

Turbulence and the k-epsilon Model Harlow and his colleagues also contributed to the early numerical modeling of turbulence and, in particular, by the postulation of the now ubiquitous k-ε model in the 60’s [51][52]. 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.

PAF (Particle-and-Force) and Free-Lagrangian Methods 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 [48]. 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

Summary of Current Research & Matrix Table

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 [49]. At any time, the particles interact with only their neighbors. The time evolution of the particles is just the solution of Newton’s equations for a multi-bodied system.

In 1965, a summary report was published [34] 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 N^2, 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 [60]. 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 70’s and 80’s in T-3 were more application-driven with close collaboration with the end users and less of explorations in CFD – an era had passed.

3. Summary of Current Research There are no obvious choices in the division of the current research in T-3. Typically a project is characterized by a unique and well-developed capability or approach, and then additional features are borrowed from other projects. Fig. 3 pictorially presents the major physics and numerical methods that are of current interest in T-3. Because of the limitation of space, the description of each project is necessarily brief. We encourage the reader to consult our Web page (http://gnarly.lanl.gov/home.html) for more information and points of contact. Where possible, readily available journals are cited as references, but when internal Los Alamos reports are cited, no outside publications were made. These reports are available through the Report Library at Los Alamos or by contacting the author.

Fig. 3. The physics and numerical approaches of current interest in T-3:

[Physics] A. Chemistry B. Turbulence C. Mixing D. Interfacial phenomena E. Phase Changes F. Solid deformation G. Electromagnetic H. Visco-elasticity I. Radiation

[Codes or Projects]

• KIVA & No-Utopia (A, B, C, D, E, 1, 3, 8, 11, 12, 13)
• Wildfire Modeling (A, B, 1, 2, 4)
• CFDLIB (A, B, C, D, G, I, 1, 2, 3, 4, 5, 6, 7, 8, 9, 12)
• FLIP (A, E, F, G, I, 1, 2, 3, 7, 8, 9)
• Granular Flow (D, F, H, 7, 9, 11)
• TELLURIDE metal casting (B, C, D, E, F, H, 1, 4, 6, 8, 13)
• Global Modeling (B, C, E, F, H, 4, 12)
• Mantle Modeling (E, F, 4, 8, 10, 12)
• Composites (D, F, H)
• Plasticity (F, H, 8, 9, 10)
• Turbulence (A, B, C, D, 3, 4, 5)

[Numerical Methods / Regimes]

1. ALE
2. Adaptivity
3. Compressible
4. Incompressible
5. Multi-field
6. Interface Tracking
7. Particles
8. Implicit
9. Lagrangian
10. Finite Elements
11. Variable Connectivity
12. Block Structured
13. Unstructured

Section 3.1 Reactive Flow and Combustion Research

3.1 Reactive Flow and Combustion Research

KIVA Family of Codes for Combustion Engines Reactive flow and combustion modeling of fully miscible species is the area of broadest use of simulation codes from T-3 and possibly from Los Alamos. Currently these are represented by the KIVA family of codes [10][11]. The KIVA codes are in worldwide use by industry, academia, and government laboratories and have a sufficient following that a yearly International KIVA Users Meeting convenes to present papers and discuss current applications and extension to the base codes. Their popularity as research tools [61], primarily because of the availability of the source code, is exemplified by the twenty-one research papers that use the KIVA codes at the 1996 SAE Spring Congress. Although the intended applications are to flow and combustion modeling in spark-ignition and diesel engines and gas turbines (as in Fig. 3.1-1), the extreme versatility and range of features have made KIVA programs attractive to a variety of non-engine applications as well. These range in scale from proposed 500-foot-high convection towers with water sprays that clean and cool the air in polluted urban areas, down to modeling silicon dioxide condensation in high pressure oxidation chambers used in the production of microchip wafers. Other applications have included the analysis of flows in automotive catalytic converters, power plant smokestack cleaning, pyrolytic treatment of biomass, design of fire suppression systems, pulsed detonation propulsion systems, stationary burners, aerosol dispersion, and design of heating, ventilation, and air conditioning systems. A complete history of KIVA as a paradigm of technology transfer from the government laboratories to industry can be found in [12].

The current version of KIVA-3 uses an unstructured mesh of hexahedrons that are groups of logical blocks of mesh and an all-speed ALE formulation from the SALE heritage. Because of the ability to model opening and closing of ports and valves, connectivity of the mesh can change during the simulation. This is a unique feature of the currently active codes in T-3 (also see CAVEAT-GT below). KIVA is also unique in that it contains a Lagrangian particle treatment of liquid spays as originally proposed by Dukowicz [36]. The current spray model includes breakup, collisions and evaporation, coupled with the turbulent gas field. This model is inherently stochastic, in contrast to the deterministic nature of all other T-3 CFD codes, and only produces an average solution for a large number of particles. The transport and chemistry equations can treat an arbitrary number of species and reactions, both kinetic and equilibrium. Mixing-controlled combustion that works in conjunction with the k-ε turbulence model and a soot model are provided.

Fig. 3.1-1. A color visualization of the temperature field during the inlet flow in a hydrogen-fueled engine simulated with KIVA-3 [62][90] . The open intake valve is on the right; the closed exhaust port is on the left. The particles are used only to simulate the flow and are colored according to the temperature of the gas. The cold intake air (blue) is shown compressing the residual hot combustion gases (red).

Parallel with the effort to continue the maturation of KIVA-3, future versions of KIVA are being developed. These codes use largely the same numerics as KIVA-3 but address the requirements of parallel computer architectures and requirements of modern mesh generation codes. KIVA-F90 is a complete rewrite of KIVA-II using Fortran 90 and executes on workstations, massively parallel architectures, and supercomputers without modification. KIVA-AC, just now under development, is an unstructured mesh version of KIVA-F90 that will support combinations of tetrahedrons and hexahedrons.

NO-UTOPIA: a Fully Implicit, Unstructured, All-Speed Flow Code This is a new reactive flow algorithm under development in T-3 for simulating a variety of hydrodynamic phenomena ranging from low-speed incompressible flows to shock-wave hydrodynamics. In-cylinder combustion of gasoline and diesel engines, external vehicle aerodynamics, underhood cooling,

Section 3.1 cont. & Section 3.2 Modern Multimaterial Methods

vehicle air conditioning, and high-velocity impact simulations are some of the applications where this algorithm may be usefully employed. This is the first three-dimensional T-3 algorithm to have a fully implicit and coupled solution to the Navier-Stokes equations, thereby eliminating the traditional separation of the Lagrangian and rezone phases. The flow solver employs an implicit finite-volume technique that performs surface integrals in parallel loops over connections associated with the two terminus nodes [79]. A continuous transition between explicit and implicit advection is employed that is second-order accurate for Courant numbers less than unity and is monotone and unconditionally stable (a 3D extension of the scheme of Collins, et al. [33]). Other unique features of the algorithm are an accurate treatment of the diffusion terms and an accurate gradient limiter for an arbitrary number of nearest neighbors associated with a node.

Wildfire Code: Crisis Response Modeling In quite a different application from reactive flows in engines, a current effort is underway in T-3 by Linn and Harlow [76] to model accurately the propagation of a forest fire. The ultimate goal is to develop a method that is sufficiently accurate but faster than real time, in order that it can be used to develop fire management strategies during a crisis. With sufficient speed, the model can be self-correcting as data is available during a crisis and can evaluate many strategies to optimize safety and effectiveness. For the method to be sufficiently fast, the description of the propagating fire must be fairly coarse, even with modern computers. Hence, the challenge is to capture accurately the effect of the expected small scale variations in the fire, such as a single burning tree, in a stochastic, global model. The techniques that are being used are the same as those used for the development of the turbulent transport equations (see the section 3.6 on Turbulence). In a similar manner, the hope is to find a universality in the scaling of the flame behavior such that the mean temperature can be related to the local burning temperatures and flame distribution. Adaptive meshing methods to high burning gradients and conformal mapping to terrain are planned features of the approach. Fig. 3.1-2 illustrates the results from the current 2D fire propagation model. This project is part of an internally funded, laboratory-wide research program to develop simulation tools for crisis management.

3.2 Modern Multimaterial Methods Even in current times, the solution of flows with large deformation of multiple materials remains a major area of research and development. In this section the tools that specialize in the interpenetration of materials are presented. Methods that treat immiscible materials are presented in the sections on modern particle methods (3.3) and the section on modeling discontinuities (3.4). The division is somewhat arbitrary, because some of the methods presented easily treat both regimes.

CFDLIB Family of Codes: Multimaterial-Multiphase Applications The numerical origins of the CFDLIB collection of codes began with a thesis study [65] in T-3 to find a stable, well-coupled, Finite Volume integration scheme for incompressible flow with co-located primitive variables. Two possibilities emerged [63][64], neither of which seemed fully satisfactory. This was because both schemes required solving two Poisson equations at each time step, rather than one Poisson solution, as in the original MAC scheme. Nevertheless, these schemes both represented cell-centered Finite Volume schemes for incompressible flow.

Fig. 3.1-2. Temperature fields of a propagating forest fire and plume 122 seconds after an intense ignition 10 m from the left boundary. A wind speed of 0.5 m/s from left-to-right propagates the fire to the right.

On another front was an internally supported project to develop a computer code, CAVEAT, for two and three dimensional, compressible flows with resolved material interfaces and large deformation that used co-located primitive variables and an ALE split computational cycle. The CAVEAT code [5] in its final incarnation included both a Godunov’s method and the so-called Total Variation Diminishing (TVD) Finite Volume method. Because of CAVEAT’s versatile block structure and highly efficient

Section 3.2 CFDLIB and Multiphase Applications

computational approach, it furnished the basic data structure that was ultimately to become CFDLIB. One of the goals of the research project was to develop a capability for integrating the compressible multiphase flow equations, and this effort focused on the TVD schemes, because Godunov’s method, an explicit scheme, had not yet been applied to multiphase flows. One of the fundamental features of TVD schemes is the use of space-time centered fluxes for advancing a cell-centered state vector. When these fluxes are exactly centered in space and time, the method is second-order accurate, and is known as the Lax-Wendroff scheme. The TVD approach was to devise a ‘limiter’ to sense when the state is tending toward new extrema, and to use the limiter to introduce a first order fluxing in such localities. During this development it became clear that classical staggered meshes and TVD space-time centered fluxes were two different ways of accomplishing the same coupling of the momentum and pressure fields. Hence, the next step was to examine a space-time centered fluxing scheme for incompressible flow. What emerged is what is now called CCMAC; a cell-centered generalization of the MAC method, which requires a single Poisson pressure solution each time step [66].

The CCMAC scheme was the key development that provided a common numerical treatment for CFDLIB: a collection of hydrocodes that are suitable for compressible flow, incompressible flow, multiphase flow of all kinds, magnetohydrodynamic forces and multi-fluid solutions, each with their own set of conservation equations. The design of each code volume in the library is modular, making the development of codes for specialized applications exceptionally fast. For example, a k-ε model of the Reynolds stress, developed for one code volume, is easily inserted into another because of the common data structure among the codes.

The FLIP approach is now being installed into CFDLIB as an option, so one can take full advantage of the nondiffusive Lagrangian approach (see section 3.3 on Particle Methods below). With the FLIP option, one can simulate the motion of a Lagrangian projectile, passing through an Eulerian gas, penetrating a Lagrangian wall, and into an Eulerian liquid.

A current area of application of CFDLIB is the modeling of a reactive flow in multiphase, multi-field problems, such as encountered in oil refining, chemicals manufacturing, metals production, and fiber processing [67]. Fig. 3.2-1 illustrates one time of a full simulation of the startup and operation of a recirculating fast-fluidized bed (FFB) reactor. Here, the goal is to model the interpenetration of gases and liquids, relative to a field of solid catalyst grains. Solid catalyst grains circulate in a flow loop consisting of a cyclone separator with a gas exit. A realistic simulation of the FFB reactor is dependent upon the physical models used to represent the effects of chemical species conversion, physical kinetics of phase change, granular flow, and multiphase fluid turbulence. The equations that embody these physical models are developed using a combination of detailed mathematics, definitive laboratory experiments, and physical intuition. The large-scale simulation is a means of bringing together these diverse sets of information in order to test the validity of the theories, and to provide important guidance to the design and operation of modern equipment.

Fig. 3.2-1. One time in a simulation of the startup and operation of a recirculating fast-fluidized bed (FFB) reactor using CFDLIB with multiple fields and phases. The colors reflect the local volume fraction of the solid catalyst (red-high, blue-low). [Labels: Vapor outlet, Cracked Oil and Catalyst, Riser Reactor, Cyclone Separator, Oil Injection Ports, Regeneration point, Catalyst]

There also exist many applications [13] that are a subset of the FFB application, such as the two-phase flow in human cardiovascular systems, or the

Section 3.3 Particle Methods: FLIP

dynamics associated with a lifeboat dropped onto the sea from a search-and-rescue aircraft [75]. These and many other contemporary applications in modern technology are addressable by the Los Alamos code library CFDLIB. Some current applications include the smelting of iron ore, alumina precipitation, combined granular and fiber flow in manufacturing, and the effects of a near-miss in the performance of defensive missiles.

3.3 Particle Methods

FLIP: a Particle Approach to Fluids The original application of the FLIP, particle-in-cell (PIC) method, was to surface ablation of laser targets. In the ablation layer, there is a transition from solid densities for which fluid-like equations are valid, to low density plasmas, for which kinetic equations are required. Particle-in-cell models of the low density plasmas had been successfully used to model the complex behavior of the low-density plasmas [43]. The goal of the FLIP development was to extend the particle-in-cell method to model fluid-like equations so that the entire ablation process could be modeled in one code.

Despite its special capabilities of PIC for following contact discontinuities, it did not give accurate solutions in general, because the transfer of information between the particles and the underlying grid resulted in numerical diffusion equal to that produced by the donor cell method. FLIP reduced this numerical diffusion to very low levels by limiting the transfer of information between grid and particles to updates [17]. That is, the particle data is updated by solutions on the grid rather than replaced. This allows the use of interpolation of arbitrary order, rather than the nearest-grid-point interpolation to which PIC was restricted. Higher order interpolation not only gives smoother solutions, but also suppresses the finite-grid-instability [20].

FLIP has several important advantages over more generally used methods. The Galilean invariance of the FLIP formulation is a distinct advantage for rotationally dominated flow problems [18], as is the conservation of angular momentum [19]. The Lagrangian description provided by the particles also can be an advantage in chemically reacting [83] and multiphase flows [92]. A mass-matrix formulation reduces numerical dissipation to extremely low levels [26]. A disadvantage of FLIP, which is apparent in adaptively zoned calculations, was eliminated by a particle density control algorithm which automatically maintains a prescribed number of particles in each cell of the mesh [71].

The principal applications of FLIP have been to modeling magnetohydrodynamic (MHD) flow in the Earth’s magnetosphere and the Sun’s heliosphere [15], to modeling plasma processing tools [72][73][74], and to modeling solids mechanics [93][94][95]. The application to plasma processing illustrates the flexibility of the method. Particles model the interaction of neutral gas and plasma with electrically active elements within a dielectric chamber as in Fig. 3.3-1. The solutions to steady-state problems are then used to calculate the motion of electrically active dust particles to aid in the design of low-contamination processing tools.

Fig. 3.3-1. Discharge in a realistic geometry using the FLIP code [74]. On the left, the electrostatic potential is shown (yellow: 50V, purple: -10V). On the right the dust density is shown, illustrating how the geometry of the base can be used to design cleaner plasma reactors.

The solid mechanics version of FLIP, uses the “material point” method, is being applied to the study of high explosives. The material point method eliminates the need for a mixture theory by evaluating equation of state and strength models particle by particle. Thus, the properties of the plastic binder, explosive granules, and gases that are evolved by the decomposition of the explosive, are each modeled accurately. Their interaction is captured by the solution of dynamical equations on a grid in the usual way.

Currently, a three-dimensional version of FLIP is being developed. A new application will be to low-gravity flow. A higher-order formulation of the immersed boundary method [92] is also under development.

Granular Flow Modeling & Section 3.4 Modeling of Discontinuities in Space

Granular Flow Modeling Aside from the application within the CFDLIB code (section 3.2), the particle approach has also been applied to one of the least understood, quite common flows, that of the flow of granular material. The knowledge of granular flow is also vital to many key industries, including chemical, material, and pharmaceutical. Lack of understanding of granular material causes serious losses in the US industry. It is estimated that eighty percent of solid processing plants encounter problems associated with granular flows.

Research on granular flow started as early as the last century, resulting in many theories focusing on the solid-like behavior of the material. Since the 1970s, the kinetic theory of molecular dynamics has been used to calculate effective stresses for large mean-free-path granular flows. In this class of theories, duration of particle collision is neglected, and only binary collisions are allowed. A systematic approach to investigate slow, dense flows has not been established. Recently, a research effort in T-3 has developed an averaged-equation method to derive the governing equations for dense flows which include interparticle contact and multiparticle interactions [101]. The advantage of the averaging method is that the closure quantities are directly computable integrals. The resulting transport equation yields, for example, a viscoelastic model for the ensemble of cohesive granular particles.

Direct numerical simulations are performed to evaluate the assumptions in the theory and to evaluate coefficients, similar to the application of direct numerical simulations for turbulence and molecular dynamics. Currently, the numerical code is used to study the shearing flow of particles coated with a thin layer of viscous liquid. With a lubrication assumption, we consider only the normal force between particles, which is modeled by a linear spring and dashpot in series. The time evolution of this system is found by solving a differential equation for each contacting particle pair [98]. For a simple shear flow, one finding is the formation of shear banding in dense, slow flows (see Fig. 3.3-2). The exact nature and physics of this is yet to be understood [77].

The averaged-equation method is a powerful method in studying granular flows and multiphase flows. An immediate generalization for the granular system is to consider the tangential force between particles, as occurs in flows of dry, granular materials. The investigation is planned of the diffusion of granular particles by the random fluctuations in the system. A key issue in understanding the unusual behavior of granular materials is the effect of microstructure and its evolution, both which can be treated with the numerical code and theory.

3.4 Modeling of Discontinuities in Space Modeling discontinuities –whether material interfaces, phase transitions, shocks or reaction fronts– have been an active area of research for T-3. Developing robust and computationally-efficient methods of numerically treating a discontinuity has seen significant progress over the decades, as in the example of shocks, from adding numerical or “artificial” viscosity to the development of approximate Riemann solutions [37] to efficiently capture shocks in realistic materials. In contrast, the robust and efficient treatment of the multidimensional aspect of discontinuities has been an outstanding challenge and has seen limited progress, despite the many applications in need. For material interfaces, the original PIC method, and later FLIP, had robustness, but not the efficiency, particularly in three dimensions, or accuracy across the interface. The “front tracking” methods have the efficiency, but suffer from extensions to three dimensions and from lack of robustness for topologically complex discontinuities. In this section, two approaches are presented for treating discontinuities, an Eulerian approach, represented by TELLURIDE and a Lagrangian approach represented by CAVEAT-GT.

Fig. 3.3-2. One time in the simulation of lubricated spheres, undergoing a shear flow from left to right, showing the formation of a shear band at an average volume fraction (58%). At the beginning of the simulation, 108 particles are placed randomly with a statistically uniform distribution in a cube with periodic boundaries.

Section 3.4 Surface Tension, VOF, and TELLURIDE

Modern Surface Tension and VOF methods Two significant advancements in modeling discontinuities in Eulerian meshes have been developed – the robust and efficient treatment of surface tension and of spatial discontinuities. The robustness of both advancements rely on the ability to neglect the complications of modeling a discontinuity below a certain spatial resolution - the dimension of the computational cell.

The continuum treatment of surface tension [21] evaluates the curvature of the phase or material discontinuity from a continuous field defined by the spatial distribution of the phases or materials. Hence, the evaluated surface tension forces are well behaved and easily calculated for three dimensional surfaces, in comparison to alternative methods.

The approach to modeling spatial discontinuities follows the tradition of the VOF method, and the current implementation is largely based on the original approach developed by Youngs [99]. The Youngs method, as applied to flows of materials, constructs a piecewise, linear (or planar in 3D) interface in each computational cell that reproduces the orientation and partial volume of a material in the cell. The reconstructed interface is then used to transport preferentially the material through the mesh. For 3D, Youngs [100] decoupled the evaluation of the normal of the interface from locating the interface within the computational cell, thereby making the three-dimensional solution tractable.

The source of the robustness of the VOF method is that the discontinuity, or interface in Youngs application, is “reconstructed” at each time increment in the simulation. By discarding the prior spatial information and constructing the discontinuity based on the new, averaged scalar information, the spatial complexity of discontinuity is limited to the mesh resolution. The contribution of T-3 to the method was to extend Youngs method of interface reconstruction on orthogonal, 2-D and 3-D meshes to 2-D, non-orthogonal and curvilinear meshes by Johnson and Zemach in the CAVEAT code [5] and later to 3-D, non-orthogonal meshes by Kothe, Zemach and Chaumeazea [69]. The implementation in CAVEAT, and later in CFDLIB, had the unique advantage that the material interface could be initially described in an accurate Lagrangian manner [81] and then could selectively be converted to a VOF treatment as the local mesh along the interface became distorted. The ability to do reconstruction in non-orthogonal cells also enables mesh adaptivity to minimize the volume of the cells that contain discontinuities, and thereby, significantly reducing the errors in the VOF approximations, which are proportional to the dimensions of the cell. A complete discussion of the history of the development of the VOF method is available [68].

TELLURIDE: A New Casting-Simulation Tool While casting-simulation tools have evolved over the past decade and are already useful in foundry operations, there are many unanswered questions, particularly those associated with the effects of the microstructural properties. Los Alamos has assembled a laboratory-wide team that is actively working on the development of a new casting-simulation computer code known as TELLURIDE. The ultimate goal is to describe virtually every property of the finished product in terms of the sequence of casting processes during its formation. The current focus is the generation of a basic numerical tool and the testing of numerous models for representing: mold-filling and fluid flow in the presence of liquid-solid phase changes (gas entrapment, surface tension effects, “slushy” flows and turbulent currents); solid-solid phase change and material response: (pure metal or alloy, grain structure, wall pull-away effects, porosity growth, bubbles or hole formation, macrosegregation, residual stress and distortion, and machinability); and geometry effects (movable walls and gate and riser design).

TELLURIDE uses Finite Volume methods with primitive variables co-located at cell centers on an unstructured hexahedral mesh and is second-order accurate in time and space. The code can execute on any modern computing platform having a Fortran-90 compiler, ranging from PCs and single workstations to symmetric multi-processor workstations, workstation clusters, and massively parallel platforms. Efficient parallelization results from explicit message passing through the use of a new parallel gather-scatter library (PGSLIB) based on the message passing interface (MPI) standard. Parallelization decreases execution time and, more importantly, allows the analyst to perform simulations with higher resolution by taking advantage of the larger memory on massively parallel platforms.

Some of the current capabilities of TELLURIDE are illustrated in Fig. 3.3-1, which display the cooling and solidification of an aluminum-copper alloy part in a mold that has already been filled and is at a uniform, initial temperature of 920°K. The underlying mesh of 5322 hexahedral cells was provided by Alcoa in a standard solids modeling format. The part is cooled from the bottom by imposing a 750°K temperature. All other surfaces are perfectly insulated. At the specified copper concentration, the cooling process

CAVEAT-GT and Section 3.5 Simulations on Parallel Platforms

includes a temperature range of 99°K over which the liquid and solid can coexist in equilibrium and the components of the alloy can migrate. Detailed information on propagation and location of the solidification front aids in the design and optimization of risers – reservoirs of excess molten alloy that seep in the part as it cools. The cooling-rate data generated by this simulation form a basis for providing insight into the final microstructure of the part. In other TELLURIDE simulations that include the injection process following the progress of the free interface using the above-mentioned surface tension and VOF techniques, the initial thermal state and alloy distribution are found to be more complex than in this idealized experiment.

Fig. 3.3-1. The liquid volume fraction (red - high, blue -low) for the cooling of a molded, alloy part. As heat is extracted from the molten part from below and solidifies, the region of solid/liquid mixture (yellow/green) moves upward.

CAVEAT-GT: a General Topology ALE method In contrast to the method presented in the last section for modeling discontinuities on Eulerian meshes, an alternative approach is to track the discontinuity, resolving the distortions of the surface to some arbitrary resolution. A development effort was undertaken in the mid 80’s to develop a general topology ALE code that could treat an arbitrary distortion of material interfaces by refining the mesh along the Lagrangian description of the interface. All prior T-3 ALE codes had been based on structured meshes, which limited the degree of distortion. It was for this capability that the VOF methods were developed for non-orthogonal meshes, as discussed in the last section. The resulting code, CAVEAT-GT [6][39], was a culmination of all the understanding gained from the earlier CAVEAT family of codes and from the Free Lagrangian efforts at Los Alamos. The numerical method was based on CAVEAT [5], but a 2D mesh of arbitrarily shaped polygons was used with piecewise-linear material interfaces. The approach included a general rezone capability that allowed the simulation to proceed in a fully Lagrangian manner until a mesh distortion criterion was exceeded and then an entirely new mesh, which preserved the material interfaces in a “near Lagrangian” fashion, was generated and then used until the next rezone was necessary. In contrast to the restrictions of a structured mesh, the general topology mesh can be chosen to have near ideal cell geometry, a regular hexagon, for minimizing numerical inaccuracies due to cell distortion. Because of the difficulty of a general rezone on a polyhedron mesh in 3D, a 3D version of the code was not contemplated.

3.5 Simulations on Parallel Platforms Throughout the history of T-3, CFD code developments have typically saturated the capabilities of existing platforms. Even so, the advent of the massively parallel platforms in the late 80’s yielded a quantum increase in speed and, even more importantly, in memory. Entire classes of problems were possible that could not even be considered earlier. And, because of the complexity of programming on these platforms, more simple CFD techniques, particularly Eulerian methods, saw a reemergence in popularity. For example, relatively simple Eulerian codes with orthogonal meshes using the VOF method described in the last section quickly displaced the more complex ALE codes, such as CAVEAT, in popularity. Now that improvements in programming tools, such as Fortran-90, and libraries of message passing routines are widely available for almost all platforms, the complexity of the CFD method is becoming less important and a greater variety of methods are being used on all platforms.

As a point of reference, earlier parallel computing strategies in T-3 focused on data parallel (SIMD) machines in which the many processors compute in lockstep. It was on these machines, such as Thinking Machines CM-200 and the CM-5, that the more simple CFD approaches excel. Later, the experience in T-3 shifted towards explicit message passing (SPMD, MIMD) machines in which the processors can execute different instructions and rely on coordinated exchange of information. This approach has the advantage of ease of transport to a variety of machines, including CM-5, nCUBE2, Cray T3D, Paragon, and, most significantly, clusters of workstations. While specialized supercomputers will continue to have important uses in the future, the area of greatest increase and availability of computing

Mantle Dynamics & Global Ocean/Climate Modeling (POP)

resources in the world will be in the area of clusters of workstations. Future CFD codes must be able to operate on these platforms.

In this section, the current code developments that have benefited from the availability of massively parallel computers are summarized. For the mantle simulations, not only did the application benefit from the new computers, but Baumgardner’s work often led the way for other applications in the development of optimal techniques and the demonstration of the importance of the new platforms.

Mantle Dynamics Simulation Thirty years after the scientific revolution in the Earth sciences that established the concept of plate tectonics, key aspects of the mechanism responsible for moving the tectonic plates across the surface of the Earth still remain unclear. Although the driving forces are generally ascribed to the process of thermal convection, the pattern of flow in the mantle responsible for the plate motions is still largely a matter of speculation. The advent of massively parallel computers, however, is now making possible high-resolution 3-D simulation of the mantle’s dynamic. When coupled with improving observational results from seismic tomography, numerical modeling of the Earth’s interior promises to yield genuine understanding of the Earth’s internal dynamics and its plate tectonics history.

The mantle dynamics code [22][23][24][25] is a 3D, finite-element model that uses a mesh comprised of hexagonal prism elements constructed from the regular icosahedron. The mantle’s silicate material is treated as a viscous fluid. A multigrid method is used to solve the global system of elliptic equations on every time step. The code utilizes domain decomposition and message-passing communication to map the calculation onto massively parallel platforms such a Cray T3D and the IBM SP-2 or onto clusters of scientific workstations. The model includes the physics of mineral-phase transition as well as temperature and stress-dependent rheology and tectonic plates over the Earth’s surface. Fig. 3.5-1 shows one time of a simulation that uses 1.3 million elements on 16 nodes of the parallel Cray T3D.

Global Ocean and Climate Modeling POP: Global Ocean Modeling The earth’s climate is determined by the complex interaction of many physical systems including the ocean, atmosphere, and biosphere. There is a growing concern that human activities (such as increasing greenhouse gases, introducing ozone-depleting chemicals, and deforestation) could alter these interactions and significantly change our current climate. However, these concerns have been very difficult to address due to our limited understanding of each of the isolated physical systems, much less the fully coupled climate system.

Fig. 3.5-1 Cutaway view of the interior of a model of the Earth’s mantle showing patterns of cold down-welling (in blue and green) and hot upwelling (in red and orange) generated by the process of thermal convection.

Performing realistic computer simulations of the global ocean is difficult because the ocean contains both relatively small spatial scales (tens of kilometers for energetic eddies) and long time scales (many centuries for the deep ocean circulation). The ocean general circulation model (OGCM), called POP (Parallel Ocean Program) [40][89], was developed to take advantage of parallel computer architectures and to perform the highest resolution global ocean simulations ever undertaken.

POP is a descendant of the Bryan-Cox OGCM which has been successful in simulating a wide range of ocean flows. The earlier model has been substantially improved and adapted for use with massively parallel computers. Improvements include a surface-pressure formulation that allows a much more realistic representation of land masses and ocean-bottom topography, an implicit free-surface technique that lets the air-sea interface evolve freely [41], and the ability to use any locally orthogonal horizontal grid

Sea Ice Model and Atmosphere/Ocean Coupling

which easily allows the Arctic ocean to be included in simulations without the problems associated with the convergence at the North Pole.

High-resolution simulations, using POP, of the global ocean on 512 nodes of the Thinking Machines CM-5 computer at Los Alamos gives a horizontal resolution from 30 km near the equator to 7 km at the poles. Estimates of the true atmospheric winds from 1985 to 1994 are used to force the model.

Sea Ice Model A component of the POP effort is a description of the sea ice (the CICE model). Many factors affect the growth and movement of sea ice, including heat fluxes from the atmosphere and ocean, short- and long-wave radiation, humidity, snowfall, albedo, salinity, ocean currents and winds [80][84]. CICE has three interacting components: a thermodynamic model that produces local growth rates of snow and ice due to heating, snowfall and local air, ice and ocean temperatures; a transport model consisting of five continuity equations for ice and snow thicknesses; and a model of ice dynamics, which predicts the velocity field of the ice pack based on a model of the material strength of the ice and forcing by wind and ocean currents. The principal difficulty in the CICE model is the description of the viscous-plastic deformation of the ice under flow. Parallel architectures present problems for these implicit models, which require excessive communication between processors. On the other hand, computationally simpler models tend to compromise important physical approximations. CICE features a new elastic-viscous-plastic model of ice dynamics that permits a fully explicit implementation and produces results equivalent to the standard implicit viscous-plastic model.

Atmosphere and Ocean Coupling A new climate model for parallel computer architectures is under development with the goal of coupling individual component models for the atmosphere, hydrosphere, sea ice and biosphere. In most existing coupled models, simulations are run at very coarse resolution (5 degrees) to enable very long integrations. The Los Alamos model combines the Community Climate Model (CCM3) from the National Center for Atmospheric Research (NCAR) for the atmosphere and biosphere with the POP ocean model in the DOE-sponsored CHAMMP (Computer Hardware, Advanced Mathematics and Model Physics) Program.

In addition, the NCAR Flux Coupler is used to “glue” the models together. The Flux Coupler accepts state variables from each component model, calculates all interface fluxes (e.g. heat and water), and sends these fluxes back to the component models. Los Alamos has added 1st and 2nd-order conservative interpolation schemes [38] that work for more general grids like the POP displaced-pole grid. Currently, the CHAMMP model is being tested with shorter simulations before beginning long (100-year) simulations to examine climate variability.

Fig. 3.5-2. This figure shows the sea surface temperature averaged over 3 years as simulated by the Parallel Ocean Program (POP).

Section 3.6 Turbulence Modeling

3.6 Turbulence Modeling: Two-Point Correlations and Engineering Models

History of Turbulence Modeling The first efforts in “turbulence” modeling directed towards a practical closure date back to Chou [28] in the early 1940’s and to Rotta [88] in the early 1950’s. These early attempts at modeling typically involved a transport equation for the turbulent kinetic energy, and met with limited success. It seemed that something was missing: a time-scale or a length-scale. One of the most successful “recipes” for producing the missing item was the introduction of the turbulent kinetic energy dissipation rate equation (“ε”), by Harlow and Nakayama in the late 1960’s [51][52]. Rather than attempt to derive the ε-transport equation, Harlow and Nakayama produced it from the turbulent kinetic energy equation and dimensional considerations alone, which lead to the still-popular, and useful, k-ε and Rij-ε [35] family of closures.

During this time period there were also attempts to develop theoretical models of turbulence, such as the Quasi-Normal models initiated by Millionshtchikov [82] in the 1940’s and the Direct Interaction Approximation of Kraichnan [70] in the 1950’s. These models differ from the “practical” closures in some profound ways. The “practical” closures to this day are based on the joint probabilities of the fluctuating fluid quantities at a single point in space and time, e.g., Rij(x,t) = ui(x,t)uj(x,t). The more fundamental theories typically consider the joint probabilities at two points, Rij(x1, x2, t) = ui(x1, t)uj(x2, t) and in some cases (e.g., DIA) at two times, Rij(x1, x2, t1, t2) = ui(x1, t1)uj(x2, t2).

Although the single point models provided a tractable set of equations, they sacrificed significant physical fidelity. But, the more fundamental closures are nearly intractable, unless severe restrictions are made (i.e., isotropy, or homogeneity), and then cast in terms of Fourier series. What do the more fundamental theories offer? The two-point models do not require a restrictive coupling of length-scales and time-scales, a type of statistical self-similarity necessary to characterize the multiscale problem by, say, two scales, k and ε. While it will probably be some time before computers are capable of solving the fundamental theories for practical problems, they provide considerable guidance in the assumptions inherent in the derivation of single-point equations, and also serve as a useful bridge between direct numerical simulations and single-point closures.

Current Research in T-3 Much of the work on turbulence in T-3 involves exploiting the more fundamental two-point turbulence modeling approach to derive “enhanced” engineering models. Thus, the effort may be considered as categorized into three subareas.

Derivation and validation of two-point models A recent development in T-3 by Besnard, Harlow, Rauenzahn and Zemach [14] has been a tractable spectral (two-point) closure by borrowing ideas from both single-point and two-point models. This model does not require an ε equation and thus does not invoke many of the self-similarity constraints implicit in the one-point closures. The model has served as the basis of a great deal of the spectral modeling work and practical simulations of turbulent mixing in multiple material problems. Extensions have been made to variable-density turbulence [32] and inhomogeneous variable-density turbulence (e.g., Rayleigh-Taylor mixing) [91].

More recently, attempts are being made in T-3, most notably by Turner, to construct an Eddy-Damped Quasi-Normal Markovian (EDQNM) model for inhomogeneous flows which does not invoke any assumptions of local homogeneity or local isotropy, and thus fundamentally differs from the work of the French researchers at L’Ecole Centrale de Lyon. The effort follows closely our direct numerical simulations using pseudo-spectral algorithms of inhomogeneous turbulence. Figure 3.6-1 shows results from the EDQNM model. Our research is directed toward understanding the degree to which the EDQNM class of models can represent the strongly intermittent zones at the edges of the turbulent zones, as well as understanding how well “gradient-diffusion” models of the one-point variety can represent the spreading of the turbulence. The work illustrates the relatively strong departures of the statistics from a near-Gaussian distribution and the importance of the action of the triple-velocity correlations on the distribution of the turbulence, and thus highlights the need, and challenge, of higher order turbulence closures.

Direct numerical simulation (DNS) and two-point closures are also being used to investigate the nature of the correlations of the fluctuating pressure-velocity and fluctuating pressure-strain. It is found that the “local” representation in differential form (rather than integral over the field) of these terms as used in most single-point closures may lead to errors at least as large as those due to truncating the hierarchy of moment equations in the “classical” closure problem. This is also true of spectral models which reduce the

Section 3.6 cont. & Section 5.7 Dynamics of Materials

vector-k-space spectral equations to a scalar k-space. By understanding the nature and extent of this “nonlocal” phenomenon, we hope to derive useable approximations that capture this feature of the pressure correlations.

Fig. 3.6-1. Results of the EDQNM model for inhomogeneous turbulence. Shown is the turbulent kinetic energy (red-high, purple-low) for the self-propagating or “diffusing” turbulence from an initial localized turbulence in the center of a channel.

Examining the two-point closures for emergent scalings and self-similarities. One of the more interesting features observed in the spectral model of Besnard et al. [14] is the emergence of self-similar spectra for turbulence undergoing homogeneous mean shears and strains and during free decay of homogeneous anisotropic turbulence at very high Reynolds numbers [29][31]. The emergence of self-similar spectra has also been observed for the case of a spectral model applied to Rayleigh-Taylor mixing by Steinkamp et al. [91]. The emergence of these self-similar spectral forms indicates that the turbulence model results (and, one hopes, the turbulence itself) can be described by far fewer degrees of freedom than required in the two-point description. The degree to which the spectral models describe actual turbulence is judged by comparisons with direct numerical simulations at low Reynolds numbers and, to whatever degree possible, by comparison to actual experiments.

One-Point “Engineering” Models The emergence of the self-similar forms suggests that the spectral models may provide a useful tool to judge the applicability of the single-point closures for a given class of flows. If the emergent self-similar form is reasonably simple, or “simplifiable”, they can be inserted into the two-point model equations, and a one-point model can be derived by the construction of appropriate integral moments. Besnard et al. have shown that if one chooses to represent the spectra as a particular self-similar form, one can then directly derive a k-ε model by constructing appropriate moments of the spectra. If the spectra from the model produce a different self-similar form in different circumstances, then a new set of moment equations may produce an improved k-ε model. An example of this has been demonstrated by Clark [29] and Clark and Zemach [30].

Current Applications The turbulence transport models developed in T-3 are being applied to examples at all flow speeds from far subsonic (incompressible) to supersonic, and with various combinations of interpenetrating fluids or clouds of droplets or grains. Some of these examples are: the fuel-air interaction in an internal combustion engine, turbulent flame behavior, unstable deformation of inertial confinement fusion capsules, nozzle flows with aerodynamic applications, two-phase flow of catalytic particles and petroleum in an industrial cracker, research problems for extended model development (free shears and mixing layers), and fluidized beds.

5.7 Dynamics of Materials In what may seem an unexpected area for a fluid dynamics group, in the past decade a significant research effort in solid mechanics has grown in T-3. The point of relevance is that many of the applications of interest involve large deformations of solids that make use of numerical methods or theoretical approaches that have been developed for fluids. The description of the dynamic response of solids to high-velocity, large deformation loadings involve strains between 0 and 200, strain rates from 10^-2 s^-1 to 10^7 s^-1, and temperatures from 300°K to beyond the melting point. Nonlinear, anisotropic inelastic material response including strain-rate phenomena, thermal softening, hardening, and failure must be considered. The resulting material models must be compatible with the incremental, continuum formulations inherent to large deformation, numerical approaches and be numerically robust and computationally efficient for large-scale computational simulations.

High Strain-Rate Plasticity The flow stress of a metal is affected strongly by the rate of deformation. Efforts are being pursued to implement plasticity models, which accurately model

Materials Dynamics, Conclusions, and Acknowledgments

hardening phenomena due to strain and strain-rate, into computational analyses [78]. Furthermore, the effects of the anisotropic plastic deformation on problems related to metal forming, machining, and impact events have been considered. Material anisotropy is a consequence of the inherent crystalline structure of a solid. Information obtained from experiments and microstructural investigations are used to construct anisotropic yield surfaces for a number of materials. The evolution of texture, including twinning phenomena, also has been pursued.

Composites There is considerable interest in using engineered composite materials to develop lighter structures that are strong under adverse conditions. Mechanical attributes, which can be tailored to provide high strength and stiffness, light weight, abrasion resistance, improved damage tolerance, and inexpensive fabrication requirements, have established composites as ideal materials for many structural applications. Composite models are under development for both epoxy and metal matrix composites and use a homogenization approach, which accounts for the response of the constituents and the interfaces within the composite [7]. These techniques provide an approach for modeling the microstructure of heterogeneous materials without requiring a detailed resolution of the micromechanics. Homogenization techniques are applicable to general loading conditions, compatible with a wide variety of constitutive models for the constituents, and are easy to implement into simulation codes [1].

Material Failure Phenomena including shear banding and the nucleation and growth of porosity and cracks must be considered for the large deformation of materials. An engineering analysis cannot provide the resolution necessary for a detailed description of these micromechanical processes. Therefore, solutions of micromechanical events are obtained in terms of the far-field distributions of the stress state. The addition of failure phenomena to a continuum analysis also can lead to numerical complexities such as strain localization. The equations governing the dynamics of failure may become ill-posed in regions of material softening. Efforts to mitigate these numerical problems have resulted in the development of rate-dependent, overstress material models for ductile failure [2][4][59] and brittle failure [3].

Conclusions The research in CFD in the early years of T-3 was prolific, even though all developments were significant at the time when so few methods were available. It was an exceptionally exciting time in the development of a new scientific field: computational fluid dynamics. The field is now fully mature, some three decades later. A comparison of early work with current research illustrates that many of the same methods are still fundamental, and modern advances in CFD are often improvements on older techniques. As a consequence, the current contributions to CFD, as represented by the current work in T-3, are directed towards developing new material models, integrating a wider variety of physics and applications to complex systems. The combination of more powerful computers with more sophisticated simulation tools have addressed problems that could not have been considered just a decade ago.

Acknowledgments The author wishes to thank the members of T-3 who both contributed much of the initial text that comprised this manuscript and who provide a truly stimulating environment and an unmatched collective resource. Special thanks to Tony Hirt who provided much material on the early days of T-3. I also wish to thank the organizers of the Canada CFD Conference for their interest and support that made this paper possible. Special thanks to M.A. Findley for the assistance in the final preparation of the manuscript.

References [1] - [26]

References [1] Addessio, F.L., and J.B. Aidun, “Analysis of Shock-Induced Damage in Fiber Reinforced Composites,” High-Pressure Shock Compression of Solids, Vol. 3, L. Davison and M. Shahinpoor, Eds., Springer-Verlag (in print). [2] Addessio, F.L., and J.N. Johnson, “Rate Dependent Ductile Failure Model,” J. Appl. Phys., Vol. 74, p. 1640, 1993. [3] Addessio, F.L., and J.N. Johnson, “A Constitutive Model for the Dynamic Response of Brittle Materials,” J. Appl. Phys., Vol. 67, p. 3275, 1990. [4] Addessio, F.L., J.N. Johnson and P.J. Maudlin, “The Effect of Void Growth on Taylor Cylinder Impact Experiments,” J. Appl. Phys., Vol. 73, p. 7288, 1993. [5] Addessio, F.L., J.R. Baumgardner, J.K. Dukowicz, N.L. Johnson, B.A. Kashiwa, R.M. Rauenzahn, and C. Zemach, “CAVEAT: A Computer Code for Fluid Dynamics Problems with Large Distortion and Internal Slip,” Los Alamos National Laboratory Report LA-10613-MS Rev. 1, May 1992. [6] Addessio, F.L., M.C. Cline, and J.K. Dukowicz, “A General Topology, Godunov Method,” Comp. Phys. Comm., Vol. 48, pp. 65–73, 1988. [7] Aidun, J.B., and F.L. Addessio, “An Enhanced Cell Model with Nonlinear Elasticity,” J. Comp. Mat., Vol. 30, p. 7248, 1996. [8] Amsden, A.A., and F.H. Harlow, J. Comp. Phys., Vol. 3, p. 80, 1968. [9] Amsden, A.A., H.M. Ruppel, and C.W. Hirt, “SALE: A Simplified ALE Computer Program for Fluid Flow at All Speeds,” Los Alamos National Laboratory Report LA-8095, 1980. [10] Amsden, A.A., “KIVA-3: A KIVA Program with Block Structured Mesh for Complex Geometries,” Los Alamos National Laboratory Report LA-12503-MS, March 1993. [11] Amsden, A.A., P.J. O’Rourke, and T.D. Butler, “KIVA-II: A Computer Program for Chemically Reactive Flows with Sprays,” Los Alamos National Laboratory Report LA-11560-MS, May 1989. [12] Amsden, D.C., and A.A. Amsden, “The KIVA Story: A Pardigm of Technology Transfer,” IEEE Trans. on Prof. Comm., Vol. 36, pp. 190–195, 1993. [13] Apen, P.G., E.L. Joyce, B.A. Kashiwa, and W.B. VanderHeyden, “Virtual National Laboratory for the Investigation of Multiphase Dynamics,” Los Alamos National Laboratory Report LA-UR-95-2968, 1996. [14] Besnard, D.C., F.H. Harlow, R. Rauenzahn and C. Zemach, “Spectral Transport Model for Turbulence,” Los Alamos National Laboratory Report LA-UR-92-1666, 1992. [15] Birn, J., F. Iinoya, J.U. Brackbill, and M. Hesse, “A Comparison of MHD Simulations of Magnetotail Dynamics,” Geophys. Res. Lett., 1995 (accepted). [16] Brackbill, J.U., and W.E. Pracht, “An Implicit, Almost-Lagrangian Algorithm for Magnetohydrodynamics,” J. Comp. Phys., Vol. 13, p. 455-482, 1973. [17] Brackbill, J.U., and H.M. Ruppel, “FLIP: A Method for Adaptively Zoned, Particle-In-Cell Calculations of Fluid Flows in Two Dimensions,” J. Comp. Phys., Vol. 65, p. 314, 1986. [18] Brackbill, J.U., “FLIP MHD: A Particle-In-Cell Method for Magnetohydrodynamics,” J. Comput. Phys., Vol. 96, pp. 163–192, 1991. [19] Brackbill, J.U., “On Modeling Angular Momentum and Vorticity in Compressible Fluid Flow,” Comput. Phys. Comm., Vol. 47, pp. 1–16, 1987. [20] Brackbill, J.U., “The Ringing Instability in Particle-In-Cell Calculations of Low-Speed Flow,” J. Comput. Phys., Vol. 75, pp. 469–492, 1988. [21] Brackbill, J.U., D. Kothe, C. Zemach, “A Continuum Method for Modeling of Surface Tension,” J. Comp. Phys., Vol. 100, pp. 335-354, 1992. [22] Bunge, H.P., and J.R. Baumgardner, “Mantle Convection and Modeling on Parallel Virtual Machines,” Computers in Physics, Vol. 9, pp. 207–215, 1995. [23] Bunge, H.P., and M.A. Richards, “The Origin of Long-Wavelength Structure in Mantle Convection: Effects of Plate Motions and Viscosity Stratification,” Geophys. Res. Lett., 1996 (submitted). [24] Bunge, H.P., M.A. Richards, and J.R. Baumgardner, “A Sensitivity Study of 3-D Spherical Mantle Convection at 10^8 Rayleigh Number: Effects of Depth-Dependent Viscosity, Heating Mode, and an Endothermic Phase Change,” J. Geophys. Res., 1996 (submitted). [25] Bunge, H.P., M.A. Richards, and J.R. Baumgardner, “The Effect of Depth-Dependent Viscosity on the Platform of Mantle Convection,” Nature, Vol. 379, pp. 436–438, 1996. [26] Burgess, B., J.U. Brackbill, and D. Sulsky, “Mass Matrix Formulation of the FLIP Particle-In-Cell Method,” J. Comput. Phys., Vol. 103, pp. 1–15, 1992.

References [27] - [54]

[27] Butler, T.D., and P.J. O’Rourke, “A Numerical Method for Two-Dimensional Unsteady Reacting Flows,” Sixteenth International Symposium on Combustion , pp. 1503–1515, 1976. [28] Chou, P.Y., “On an Extension of Reynolds’ Method of Finding Apparent Stress and the Nature of Turbulence,” Chin. J. Phys., Vol. 4, pp. 1–33, 1940. [29] Clark, T.T., “Spectral Self-Similarity of Homogeneous Anisotropic Turbulence,” Ph.D. Dissertation, University of New Mexico, 1991, also available as Los Alamos National Laboratory Report LA-12284-T, 1992. [30] Clark, T.T., and C. Zemach, “Engineering Turbulence Models from the Standpoint of Spectral Turbulence Models,” Bull. Am. Phys. Soc., Vol. 37 (8), p. 1799, 1992. [31] Clark , T.T., and C. Zemach, “A Spectral Model Applied to Homogeneous Turbulence,” Phys. Fluids, Vol. 7 (7), pp. 1674–1694, 1995. [32] Clark, T.T. , and P.B. Spitz, “Two-Point Correlation Equations for Variable-Density Turbulence,” Los Alamos National Laboratory Report LA-12672-MS, 1995. [33] Collins, J. P., P. Colella, H.M. Glaz, “An Implicit-Explicit Eulerian Godunov Scheme for Compressible Flows,” J. Comp. Phys., Vol. 116, pp. 195-211, 1995. [34] Daly, B.J., F.H. Harlow, and J.E. Welch, “Numerical Fluid Dynamics Using the Particle-And-Force Method,” Los Alamos National Laboratory Report LA-3144, 1965. [35] Daly, B.J., and F.H. Harlow, “Transport Equations in Turbulence,” Phys. Fluids, Vol. 13, pp. 2634–2649, 1970. [36] Dukowicz, J.K., “A Particle-Fluid Numerical Model for Liquid Sprays,” J. Comp. Phys., Vol. 35, pp. 229-253, 1980. [37] Dukowicz, J.K., “A General, Non-Iterative Riemann Solver for Godunov’s Method,” J. Comp. Phys., Vol. 61, pp. 119–137, 1985 [40]. [38] Dukowicz, J.K., and J. Kodis, SIAM J. Sci. Stat. Comp. , Vol. 8, pp. 305-321, 1987. [39] Dukowicz, J.K., M.C. Cline, and F.L. Addessio, “A General Topology Godunov Method,” J. Comp. Phys., Vol. 82, pp. 65–73, 1989. [40] Dukowicz, J.K., R.D. Smith, and R.C. Malone, “A Reformulation and Implementation of the Bryan-Cox-Semtner Ocean Model on the Connection Machine,” J. Atmos. Oceanic Technol., Vol. 10, pp. 195–208, 1993. [41] Dukowicz, J.K., and R.D. Smith, “Implicit Free-Surface Method for the Bryan-Cox-Semtner Ocean Model,” J. Geo. Res., Vol. 99, pp. 7991-8014, 1994. [42] Evans, M.W., and F.H. Harlow, “The Particle-In-Cell Method for Hydrodynamic Calculations,” Los Alamos National Laboratory Report LA-2139, 1957. [43] Forslund, D.W., and J.U. Brackbill, “Magnetic Field Induced Surface Transport in Laser Irradiated Foils,” Phys. Rev. Let., Vol. 48, p. 1614, 1982. [44] Fromm, J. and F.H. Harlow, “Numerical Solution of the problem of Vortex Sheet Development,” Phys. Fluids, Vol. 6, p. 975, 1963. [45] Gentry, R., R.E. Martin and B.J. Daly, J. Comp. Phys., “An Eulerian Differencing Method for Unsteady Compressible Flow Problems,” Vol. 1, p. 87, 1966 [46] Harlow, F. H., “Hydrodynamic Problems Involving Large Fluid Distortion,” J. Assoc. Comp. Mach., Vol. 4, p. 137, 1957. [47] Harlow, F.H., D.O. Dickman, D.E. Harris, and R.E. Martin, “Two-Dimensional Hydrodynamic Calculations,” Los Alamos National Laboratory Report LA-2301, 1959. [48] Harlow, F.H., and B.D. Meixner, “The Particle-And-Force Method for Fluid Dynamics,” Los Alamos National Laboratory Report LA-MS-2567, 1961. [49] Harlow, F.H., “Theory of Correspondence between Fluid Dynamics and Particle-And-Force Method,” Los Alamos National Laboratory Report LAMS-2806, 1963. [50] Harlow, F.H., and E. Welch, “Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluids with Free Surface,” Phys. Fluids, Vol. 8, p. 2182, 1965. [51] Harlow, F.H., and P.I. Nakayama, “Turbulence Transport Equations,” Phys. Fluids, Vol. 10, pp. 2323–2332, 1967. [52] Harlow, F.H., and P.I. Nakayama, “Transport of Turbulence Energy Decay Rate,” Los Alamos National Laboratory Report LA-3854, 1968. [53] Harlow, F.H., A.A. Amsden, and C.W. Hirt, 2nd Int. Conf. on Num. Methods of Fluid Dynamics, Berkeley, CA, Sept. 1970. [54] Hirt, C.W., A.A. Amsden, and J.L. Cook, “An Arbitrary Lagrangian-Eulerian Computing Method for All Flow Speeds,” J. Comp. Phys., Vol. 14, p. 227, 1974.

References [55] - [80]

[55] Hirt, C.W., B.D. Nichols, and N.C. Romero, “SOLA—A Numerical Solution Algorithm for Transient Fluid Flows,” Los Alamos National Laboratory Report LA-5852, 1975. [56] Hirt, C.W., N.C. Romero, M.D. Torrey, and J.R. Travis, “SOLA-DF: A Solution Algorithm for Nonequilibrium Two-Phase Flow,” Los Alamos National Laboratory Report LA-7725-MS (and NUREG/CR-0690), 1979. [57] Hirt, C.W., and B.D. Nichols, “Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries,” J. Comp. Phys., pp. 39, 201, 1981. [58] Hirt, C.W., J. Comp. Phys., Vol. 2, p. 339, 1968. [59] Johnson, J.N., and F.L. Addessio, “Tensile Plasticity and Ductile Fracture,” J. Appl. Phys., Vol. 64, p. 6699, 1988. [60] Johnson, N.L., “Discrete Macroscopic Simulation of Fluid Flows,” Doctoral Thesis, University of Wisconsin, Madison, 1983. [61] Johnson, N.L., and Y. Takagi, editors. Modeling in Diesel and SI Engines, SAE publication SP-1123, 1995. [62] Johnson, N. L., A. A. Amsden, T. D. Butler, “Hydrogen Program Combustion Research: Three Dimensional Computer Modeling,” Proceedings of the 1995 DOE/NREL Hydrogen Program Review, Vol. II, pp. 285-304, NREL/CP-430-20036, 1995. [63] Kashiwa, B.A., “A Generalized MAC Method for Incompressible Fluid Flow,” Los Alamos National Laboratory Report LA-10853-MS, November 1986. [64] Kashiwa, B.A., “A Predictor-Corrector Method for Calculation of Incompressible Fluid Flow Dynamics,” Los Alamos National Laboratory Report LA-10854-MS, November 1986. [65] Kashiwa, B.A., “Statistical Theory of Turbulent Incompressible Multimaterial Flow,” University of Washington Doctoral Thesis, Los Alamos National Laboratory Report LA-11088-T, 1987. [66] Kashiwa, B.A., N.T. Padial, R.M. Rauenzahn, and W.B. VanderHeyden, “A Cell-Centered ICE Method of Multiphase Flow Simulations,” in Proceedings, ASME Symposium on Numerical Methods for Multiphase Flows, Lake Tahoe, NV, 19–23 June 1994. [67] Kashiwa, B.A., “Reactive Multiphase Flow Simulation Workshop Summary,” W.B. VanderHeyden, editor, Los Alamos National Laboratory Report LA-UR-95-2194, June 1996. [68] Kothe, D., and W. Rider, “Comments on Modeling Interfacial Flows with Volume-of-Fluid Methods,” Los Alamos National Laboratory, LA-UR-94-3384, 1994, submitted to J. Comp. Phys. , Feb. 1995. [69] Kothe, D., C. Zemach and A. Chaumeazea, manuscript in preparation, 1996. [70] R. H. Kraichnan, “The Structure of Isotropic Turbulence at Very High Reynolds Number,” J. Fluid Mech. , Vol. 5, pp. 497-543, 1959. [71] Lapenta, G., and J.U. Brackbill, “Control of the Number of Particles in Fluid and MHD Particle-In-Cell Methods,” Comp. Phys. Comm., Vol. 87, p. 139, 1995. [72] Lapenta, G., F. Iinoya, and J.U. Brackbill, “Particle-In-Cell Simulation of Glow Discharges in Complex Geometries,” IEEE Transactions on Plasma Science, Vol. 23, p. 769, 1995. [73] Lapenta, G., and J.U. Brackbill, “Immersed Boundary Method for Plasma Simulation in Complex Geometries,” IEEE Trans. Plasma Sci., 1996 (accepted). [74] Lapenta, G., and J.U. Brackbill, “Simulation of Dust Particle Dynamics for Electrode Design in Plasma Discharges,” Plasma Sources Science and Technology, 1996 (submitted). [75] Lewis, M.W., B.A. Kashiwa, R.W. Meier, and S. Bishop, “Nonlinear Dynamic Fluid-Structure Interaction Calculations with Coupled Finite Element and Finite Volume Programs,” ASME Winter Annual Meeting, Chicago, IL, 1994. [76] Linn, R.R., and F.H. Harlow, “Wildfire Propagation,” Los Alamos National Laboratory Report, LALP-96-22, 1996. [77] Lun, C.K.K., and A.A. Bent, “Numerical Simulation of Inelastic Frictional Spheres in Simple Shear Flow,” J. Fluid Mech., Vol. 258, pp. 335–353, 1994. [78] Maudlin, P. J., S.I. Wright, G. T. Gray III, and J. W. House, “Application of Faceted Yield Surfaces for Simulating Compression Tests of Textured Materials,” Metallurgical and Materials Applications of Shock-Wave and High Strain Phenomena, editors L. E. Murr, K.P. Staudhammer and M.A. Meyers, Elsevier Science Publishers, 1995. [79] Mavriplis, D.J., “Three-Dimensional Multigrid Reynolds-Averaged Navier-Stokes Solver for Unstructured Grids,” AIAA Journal, Vol. 33, pp. 445-453, 1995. [80] Maykut, G.A., and N. Untersteiner, “Some Results from a Time Dependent Thermodynamic Model of Sea Ice,” J. Geophys. Res., Vol. 76, pp. 1550–1575, 1971.

References [81] - [101]

[81] Menikoff, R., K. Lackner, N.L. Johnson, S.A. Colgate, J.M. Hyman, G.A. Miranda, “Shock Driven by a Phased Implosion,” Phys. Fluids A, Vol. 3, pp. 201–218, 1991. [82] Millionshtchikov, M., “On the Theory of Homogeneous Turbulence,” Dokl. Akad. Nauk. SSSR, Vol. 32, pp. 615–618, 1941. [83] O’Rourke, P.J., B. Larratirou, and J.U. Brackbill, “On Particle-Grid Interpolation in Particle-In-Cell Methods and Application to a Model Flame Problem,” J. Comp. Phys., Vol. 109, pp. 37–52, 1993. [84] Pounder, E.R., “The Physics of Sea Ice,” in The Sea, M.N. Hill, Ed., Vol. 1, pp. 826–838, Interscience, 1962. [85] Ramshaw, J.D., and J.K. Dukowicz, “APACHE: A Generalized-Mesh Eulerian Computer Code for Multicomponent Chemically Reactive Fluid Flow,” Los Alamos National Laboratory Report LA-7427, January 1979. [86] Rivard, W.C., and M.D. Torrey, “K-FIX: A Computer Program for Transient, Two-Dimensional, Two-Fluid Flow,” Los Alamos National Laboratory Report LA-NUREG-6623, 1976. [87] Rivard, W.C., O.A. Farmer, and T.D. Butler, “RICE: A Computer Program for Multicomponent Chemically Reactive Flows at All Speeds,” Los Alamos National Laboratory Report LA-5812, 1975. [88] Rotta, W., “Statische Theorie Nichthomogener Turbulenz,” Z. Fluids, Vol. 129, p. 547, 1951. [89] Smith, R. D., J. K. Dukowicz, and R. C. Malone, “Parallel Ocean General Circulation Modeling,” Physica D Amsterdam, Vol. 60, pp. 38-61, 1992. [90] J. R. Smith, S. Aceves, N. L. Johnson, and A. A. Amsden. “Progress towards an Optimized Hydrogen Series Hybrid Engine.” Alternative Fuels and Natural Gas, ASME ICE-Vol. 25-3; edited by J. A. Caton, 1995. [91] Steinkamp, M.J., T.T. Clark, and F.H. Harlow, “Stochastic Interpenetration of Fluids,” Los Alamos National Laboratory Report LA-13016-MS, 1995. [92] Sulsky, D., and J.U. Brackbill, “A Numerical Method for Suspension Flow,” J. Comp. Phys., Vol. 96, pp. 339–368, 1991. [93] Sulsky, D., Z. Chen, and H.L. Schreyer, “The Application of a Material-Spatial Numerical Method to Penetration,” in New Methods in Transient Analysis, ASME, pp. 91–102, 1992. [94] Sulsky, D., and H.L. Schreyer, “A Particle Method with Large Rotations Applied to the Penetration of History-Dependent Materials,” in Symposium on Advances in Numerical Simulation Techniques for Penetration and Perforation of Solids, The American Society of Mechanical Engineers, 1993. [95] Sulsky, D., and H.L. Schreyer, “The Particle-In-Cell Method as a Natural Impact Algorithm,” Sandia National Laboratories, Contract No. AC-1801, 1993. [96] Travis, J.R., C.W. Hirt, and W.C. Rivard, “Multidimensional Effects in Critical Two-Phase Flow,” Nuc. Sci. Eng., Vol. 68, p. 338, 1978. [97] Viecelli, A., J. Comp. Phys., Vol. 4, p. 543, 1969. [98] Walton, O.R., and R.L. Braun, “Stress Calculations for Assemblies of Inelastic Spheres in Uniform Shear,” Acta Mech., Vol. 63, pp. 73-86, 1986. [99] Youngs, D. L., “Time-Dependent Multi-Material Flow with Large Fluid Distortion,” in Numerical Methods for Fluid Dynamics, edited by K. W. Morton and M. J. Baines, pp. 273-285, 1982. [100] Youngs, D. L. , “An Interface Tracking Method for a 3D Eulerian Hydrodynamics Code,” Technical Report 44/92/35, AWRE, 1984. [101] Zhang, D.Z., and A. Prosperetti, “Averaged Equations for Inviscid Disperse Two-Phase Flow,” J. Fluid Mech., Vol. 267, pp. 185-219, 1994.