Journal of Computational Physics 195 (2004) 413 www.elsevier.com/locate/jcp Short Note Introduction to HarlowÕs scientiﬁc memoir J.U. Brackbill T-3, Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA Received 3 September 2003; accepted 3 September 2003 In physics and mathematics, the names we give important contributions remind us of their origin, such as Schr€ odingerÕs equation or the Lax equivalence theorem. In scientiﬁc computation, this custom is sometimes followed, as in the Metropolis algorithm or GodunovÕs method, and sometimes not. For example, the names of the MAC (Marker-and-Cell) and PIC (Particle-in-Cell) methods sound more like the product of a corporation than an invention of an individual. Nevertheless, these widely used methods, and many others besides, were contributed by Francis Harlow, as described in the Special Paper, ÔFluid Dynamics in T-3. . .Õ, which follows. As is Sergei GodunovÕs earlier paper , FrankÕs memoir is a ﬁrst-person account and personal in tone. It successfully conveys a sense of a time when scientiﬁc computing was more a dream than a reality, and opportunities for new discoveries lay in every direction. Frank has chosen to recount the activities of T-3, a group he formed in the Theoretical Division of Los Alamos National Laboratory, as it used large-scale computation to solve fundamental ﬂuid dynamics problems. The groupÕs publications did much to introduce the wider scientiﬁc community to the potential of computational ﬂuid dynamics. Calculations such as free-surface ﬂow resulting from a dam breaking, and of a Karmann vortex street attracted a wide audience, and laboratory reports taught how to do it themselves. FrankÕs paper also conveys the scope of T-3Õs contributions, which span compressible and incompressible ﬂow, free-surface and multiﬂuid ﬂow, and turbulence. To all of these, Frank and his group made pioneering contributions. Applied to Frank, unusual is an understatement. He is singular. I have known him for almost 40 of the 50 years he has been a staﬀ member at Los Alamos, so mine is not an entirely objective point of view. However, even a completely detached observer would be struck by the breadth of his interests. His paintings line the walls of the Laboratory, his books on Southwest Indian pottery ﬁll a shelf, and his monograph on brachiopods, in which he identiﬁes 24 new species and 3 new genera, is still the deﬁnitive reference on the lowest parts of the Pennsylvanian period. He has grown orchids, learned silversmithing, and stopped riding a Harley because it was not fast enough. He knows an incredible number of people in all walks of life by their ﬁrst name, yet can be met in person only in Los Alamos because he refuses to travel. It all seems to work for him, because he brings to his work an undiminished energy and focus. Reference  S.K. Godunov, Reminiscences about diﬀerence schemes, J. Comput. Phys. 153 (1999) 6–25. E-mail address: [email protected]v (J.U. Brackbill). 0021-9991/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2003.09.030 Journal of Computational Physics 195 (2004) 414–433 www.elsevier.com/locate/jcp Review Fluid dynamics in Group T-3 Los Alamos National Laboratory q (LA-UR-03-3852) Francis H. Harlow * Los Alamos National Laboratory, MS-B216, Los Alamos, NM 87544, USA Received 12 June 2003; received in revised form 25 August 2003; accepted 3 September 2003 Abstract The development of computer ﬂuid dynamics has been closely associated with the evolution of large high-speed computers. At ﬁrst the principal incentive was to produce numerical techniques for solving problems related to national defense. Soon, however, it was recognized that numerous additional scientiﬁc and engineering applications could be accomplished by means of modiﬁed techniques that extended considerably the capabilities of the early procedures. This paper describes some of this work at The Los Alamos National Laboratory, where many types of problems were solved for the ﬁrst time with the newly emerging sequence of numerical capabilities. The discussions focus principally on those with which the author has been directly involved. Ó 2003 Elsevier Inc. All rights reserved. Keywords: Computational ﬂuid dynamics; History of computing; Turbulence; Multi-ﬁeld ﬂow; Strong distortions; Incompressible ﬂow; Relativistic ﬂuids 1. Introduction The year 1953 saw the arrival of the ﬁrst big commercial computers at the Los Alamos Scientiﬁc Laboratory (now the Los Alamos National Laboratory). My employment in the Theoretical Division during that same year gave to a group of us the opportunity to become pioneers in the largely unexplored subject of numerical ﬂuid dynamics. Several years before my arrival, the MANIAC computer had been developed in our Division, paving the way for the innovative creation of important new features of both hardware and software during the initial stages of the computer revolution, but I was not a part of that project. q Preparation of this paper was performed under the auspices of the National Nuclear Security Administration of the US Department of Energy at Los Alamos National Laboratory under Contract No. W-7405-ENG-36. * Tel.: +1-505-667-9090; fax: +1-505-665-5926. E-mail address: [email protected] (F.H. Harlow). 0021-9991/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2003.09.031 F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 415 The ﬁrst of the large computers to arrive here was the Model 701 from the International Business Machines (IBM) Corporation. This machine was great and marvelous for all of us and we were very excited to follow all sorts of exploratory pathways for the development of new techniques to solve hard problems. The 701 was a big monster; it occupied a huge room ﬁlled with cathode-ray tubes that contained the memory of 2048 36-bit words or 4096 16-bit halfwords. It was cumbersome, and it had troubles, but we loved it. We would come at any time in the night or day to run it ourselves. It was a hands-on operation. It was something that was very special to be involved with. We did not realize what a dawn of the computing era we were witnessing. At that time we programmed in machine language without the advantages of index registers. Both the address arithmetic and the problem-solving calculations were performed by ﬁxed-point manipulations, although sometimes we used a subroutine called the Dual Interpretive Routine, which was very slow and cumbersome, but we could invoke it at selected places in our programs for doing ﬂoatingpoint arithmetic. By the time that Fortran arrived, we were so accustomed to the machine-language programming techniques that we felt suspicious of the new-fangled ‘‘black-box’’ approach and kept to the familiar tried-and-true procedures until well into the 1960s. Soon after my arrival at Los Alamos it became apparent that ﬂuid dynamics was the ﬁeld for me. My main specialty at ﬁrst was to look at what seemed very natural to me, namely supersonic ﬂows, in which it is possible to imagine step-by-step the progress of sound waves across a computational mesh, and of the shocks that go through the ﬂuid around various kinds of obstacles. I cut my teeth on supersonic ﬂows with incompressible ﬂows coming later. There is much that could be told about the pioneering work in our Group, but I have limited the discussion to include only the topics with which I was directly involved, in collaboration with many creative and talented associates as described in the Acknowledgment section of this paper, and listed as co-authors in my publications in physics and mathematics [1–151]. One of the things about supersonic ﬂows that I soon learned was that there were at that time two main ways to think about zoning space for resolving the behavior of a ﬂuid. One of them we call the Lagrangian approach, which means having a mesh of computational cells that follow the motion of a ﬂuid through its contortions and whatever processes that take place. The mesh could follow interfaces beautifully. The other, the Eulerian viewpoint, has a ﬁxed mesh of cells that stay in one place in the laboratory frame and the ﬂuid ﬂows through it. This second approach is great for looking at large distortions; there is no movable computational mesh to get tangled. But on the other hand, it has problems if you want to follow a sharp interface between two ﬂuids. Eulerian techniques tend to diﬀuse the interface and to smear out sharp shocks. (In those days we had only begun to experiment with very primitive interface reconstruction techniques.) 2. Particles in cells (1955) These problems with each of the computational viewpoints led us early in our work to think about some way to combine the two. Martha Evans and I started spending part of our time exploring ideas for accomplishing this amalgamation, and we developed a technique that turned out, despite a lot of initial skepticism, to be very useful. This technique we call the Particle-in-Cell (PIC) method [3,4,9,16,19,26]. It is a computational technique in which there is an Eulerian set of cells together with Lagrangian marker particles that move through the cells. Each marker particle has a constant mass. At that stage of the development the particles did not carry velocity or energy; when a particle moved across a cell boundary its mass was subtracted from the donor cell and added to the receptor cell, and the event also signaled that momentum and energy should be transferred. For each of the transported quantities the amount was the mass weighted fraction of the totals in the donor cell. This process can be viewed in two distinctly diﬀerent ways. One of them is based on the idea of time splitting, in which a computational cycle is broken into two parts, the ﬁrst being an update of the quantities 416 F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 by all the processes except advection, while the second part moves the particles and accomplishes all of the advective ﬂuxing. The second way, as recognized by Eleazer Bromberg (a visitor to our Group from the Courant Institute at New York University), is to consider that the mesh of cells is actually Lagrangian, so that the ﬁrst part of the calculations in a cycle advances all the variables, including the mesh coordinates and the particle positions, and the second part maps the mesh back into its original conﬁguration, leaving the particles at their new locations . Some years later, Tony Hirt and Dan Butler of our Group discovered that instead of a complete mapping back to the original mesh each cycle, the mapping could be partial, leading to a technique that is called an Arbitrary–Lagrangian–Eulerian (ALE) method, which limits to either of the two viewpoints. An interesting consequence of the transport of momentum and energy in mass-weighted proportions of the amounts in the donor cell is that the process introduces an eﬀective artiﬁcial viscosity that serves to remove some of the numerical instabilities that otherwise would have to be removed by means of an explicit artiﬁcial viscosity. The strength of the eﬀective viscosity was calculated by means of a truncation-error analysis. The double-mesh zoning, the time splitting ideas, the donor-cell concept, and the use of truncation-error analysis all subsequently developed into procedures that have received wide usage in the world of computer ﬂuid dynamics. The early versions of the PIC method had various restrictions and diﬃculties in its application to the problems that we wanted to solve. Principal among these was the necessity to use a fairly large number of particles in each cell in order to reduce the ﬂuctuations that occur as a result of particle passage across a boundary. Cells that have completely emptied create troubles for their neighbors, so that it was important to avoid that circumstance if possible, or else invent some type of remediation by which to avoid the interpretation that an empty cell represents a region of vacuum. Bart Daly, Dan Butler, and Tony Amsden made major contributions towards the remediation of these diﬃculties. In 1966 Tony Amsden wrote a report that is still a major resource for the early PIC techniques . We did not have the computer memory and speed that are available today, so that at ﬁrst we solved only relatively simple problems in two space dimensions. We also discovered that with the earliest formulation, it was necessary to restrict the time increment per cycle to half of what is allowed by the velocity-limiting part of the advective stability constraint, in order to prevent particles from crossing over each other. Because the formulation was completely explicit, we also were limited to consider only the problems that had no regions with low Mach number other than those that had been initialized to uniform pressure and zero velocity. Cells having more than one material were diﬃcult to treat correctly if each material had a complicated equation of state, so that all of our early test calculations were performed using polytropic-gas descriptions allowing us to simply sum the partial pressures. Conservation of linear momentum was ensured by the way in which the equations were formulated, but conservation of energy was an issue that revolved around whether the equations should transport internal energy or total energy. Many of these issues have subsequently been resolved. To demonstrate the validity of the PIC method for solving complex problems we investigated a variety of circumstances involving the strongly distorting dynamics of several materials driven by the passage of a shock over the conﬁguration . Examples include: passage of a shock over a bubble, the compression of rectangular slabs of one material imbedded in another material, the refraction of a shock over an oblique interface (for which there is good experimental data with which to compare), and the instability of a perturbed interface between two materials. Additional examples are described in [5,7,8,15,26]. 3. Low-speed ﬂows The PIC method was developed for application to ﬂuid ﬂows that are subjected to large distortions involving signiﬁcant compression during the course of the problem. Designed for ﬂows that have Mach F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 417 number greater than about 0.1, the ﬁnite-diﬀerence equations could be solved explicitly. For calculations of ﬂow at lower Mach numbers, it was at ﬁrst speculated that good results could be obtained from the compressible-ﬂuid codes by simply decreasing the time step per cycle. This procedure is not only more expensive in computer time, but also has signiﬁcant diﬃculties in obtaining accuracy as the ﬂow approaches complete incompressibility. For this type of problem we developed two very diﬀerent methodologies. One of these techniques is based on the properties of the vorticity and stream function. Jacob Fromm in our Group investigated the idea that the evolution of these quantities could be calculated numerically, and he soon worked out a procedure that was highly successful. The approach works especially well for conﬁned ﬂows in two space dimensions, where the vorticity vector always points in the third direction, and the boundary conditions can be expressed as constraints on the velocity. The transport for vorticity is solved explicitly for each time cycle of the calculation, and the negative of the updated value is used as the source term in the Poisson equation for the stream function. In the early days of this technique we used a Jacobi or Gauss–Seidel relaxation procedure to solve the equation. Each cycle ended with the extraction of the new velocities from the stream function. The most interesting of the output for a calculation was the pattern of time varying streamlines, obtained from contours of the stream function. Our most famous application of this technique was to the investigation of the von Karman vortex street generated by the ﬂow of an incompressible ﬂuid by an obstacle, which in those days was constrained to be rectangular. For low Reynolds numbers we obtained a steady ﬂow pattern with two downstream vortices, but at higher ﬂow speeds the street emerged and had the correct shedding frequency. We were ecstatic, and published these ﬁrst-in-the-world calculated results, which were subsequently reprinted in a Japanese reprint series [18,22,23]. The other of the two techniques that we developed for incompressible ﬂow works directly with the primitive variables of velocity and pressure. The principal incentive for this alternative development was to have a procedure that allows for the use of pressure boundary conditions, for applicability to free-surface ﬂows in both two and three dimensions. We called the resulting method the Marker-and-Cell (MAC) technique . It uses marker particles to show the ﬂuid contortions and to follow the changing position of the free surface, at which we could apply pressure boundary conditions (later generalized to include the full stress conditions). Eddie Welch and I realized that the principal challenge for developing a numerical technique to solve incompressible ﬂuid ﬂows using the primitive variables is to ﬁnd a way to ensure that the divergence of the velocity ﬁeld vanishes. For the ﬁrst version of the MAC method we attempted to describe the ﬁelds of velocity and pressure with co-located (cell-centered) values in a mesh of rectangular cells. The procedure was to be completely Eulerian, with the ﬂuid ﬂowing through the stationary mesh. We knew that a crucial requirement was to achieve the incompressibility constraint through the determination of pressures in each computational cycle consistent with the desired objective. In eﬀect the achievement of that goal also will remove the sound-speed restrictions on the time step per cycle; sound signals could propagate to all parts of the mesh instantaneously. The only problem with the idea soon turned out to be the diﬃculty in expressing the divergence constraint with cell-centered placement of velocities. To overcome this impediment we changed to placement of velocity components at the edges of each computational cell. In our ﬁrst two-dimensional version, the four sides of each cell received the normal component of velocity, placed in the center of the edge. This arrangement has subsequently come to be known as a staggered mesh. To derive the equation for pressure we used the well-established trick of taking the divergence of the moment equation, and thus obtained a Poisson equation for the pressure, which in those days we solved by using the Jacobi method to relax the ﬁeld. We also tried a Gauss–Seidel relaxation procedure, ﬁnding faster convergence, but at the price of altering the vorticity, which, during the iteration, is conserved to within round-oﬀ with the slightly slower Jacobi technique. Later, Tony Hirt discovered that a simple Newton–Raphson technique works much better and faster, at the same time improving considerably the process of inserting the desired boundary conditions on pressure and velocity. 418 F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 In eﬀect the MAC-method calculations take place in a time splitting process, with one part consisting of updating the velocities in response to all of the terms (body forces, advection, viscosity, etc.) except the pressure gradients, and the second part accomplishing the calculation of pressures that return the velocity to a state with vanishing divergence but without changing the vorticity that was updated in the ﬁrst part. After some initial struggling, we found that the MAC method worked very well, and we tried solving a wide variety of problems with results that exceeded our most optimistic expectations. One example that we calculated was the same von Karman vortex street that we had studied with the vorticity-and-streamfunction method. An interesting discovery was that the advection of velocity should not be accomplished with the donor-cell method that we often used in order to give greater numerical stability to the calculations; that procedure simply produced too much eﬀective artiﬁcial viscosity so that the desired instability would never occur. To overcome this diﬃculty, we used instead the second-order centered advection, knowing that one of the adverse eﬀects is to allow for the growth of a slowly developing numerical instability. It turned out that the instability was not overly detrimental to the calculation except in the region upstream of the obstacle. It was discovered only later by Tony Hirt that the presence of velocity gradients has a signiﬁcant eﬀect on the nature of the second-order advective scheme that we used, reducing the instability in regions of deceleration and enhancing it where the ﬂow is increasing in speed. With this new MAC technique we explored a large variety of ﬂow problems: the broken-dam conﬁguration in which the instantaneous removal of the dam frees the water to rush forward and interact with an obstacle with such force as to splash over the top; the opening of a sluice gate so that water could rush out with such power as to create a backward-breaking wave; the interaction of a jet of water impinging into a pool of water on a ﬂat plate; the formation of a hydraulic jump, or bore; the generation of a waterfall; the rise of a bubble in a ﬂuid; Rayleigh–Taylor instability carried to large amplitudes with the development of the classical ‘‘mushroom’’ conﬁgurations at the tip of the downward moving spike; the splash of a liquid drop into a pool of ﬂuid. These results were being obtained for the ﬁrst time in the world, and we were pretty excited to see them and to publish them in both Laboratory reports and refereed journals [28– 30,33,35,37,49,50]. 4. All-speed ﬂows (1971) The next challenge that we faced was to combine the explicitness of the high-ﬂow-speed techniques with the implicitness of the low-ﬂow-speed techniques, so as to obtain a general procedure that could work well for ﬂows with both magnitudes of speed. Tony Amsden and I realized that the key to the development lay in specifying the quantity that we wanted to bring to zero in a generalized Newton–Raphson iteration scheme. As it turned out there were two appropriate quantities that satisﬁed our needs. At ﬁrst we used the equation for conservation of mass, including any possible sources and sinks. Later, Tony Hirt used the diﬀerence between the pressure and the equation-of-state expression for pressure as a function of density and speciﬁc internal energy. The essence of a computational cycle using the ﬁrst convergence criterion consisted of a purely explicit update of velocities and internal energies and an iterative solution for determining the modiﬁed new densities and velocities that satisﬁed the equation for conservation of mass. The ﬁrst criterion was incorporated by Tony Amsden into a code called SIERRA and tested on some problems in which the overall ﬂow was far subsonic but the eﬀects of sound-signal propagation was essential to represent. We published the method and results, but soon realized that there were signiﬁcant improvements that could be made and accordingly wrote a second paper to describe these modiﬁcations [40,48,53]. (If anyone thinks that all of our developments took place in a direct fashion with no glitches along the way, I have to assure you that lots of false starts and turns occurred in virtually every project that we undertook.) Use of the second criterion for convergence worked well also, as long as the ﬂow was not too far subsonic. F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 419 This new method for combining all ﬂow speeds into one technique is called the ICE method. It has been exciting to see the numerous extensions that these ideas have fomented. Applications have been accomplished with a variety of new codes developed in our Group and elsewhere in the world for problems involving internal combustion engines, multiphase ﬂows, national defense issues, nuclear reactor behavior, climate, etc. Both Eulerian and Lagrangian versions have been successfully used, and other exotic numerical techniques have also incorporated these ideas. 5. Publications Many of these results were published in reports and journal articles (see Bibliography), but not without some opposition. There was, for example, a conversation that was reported to me, which took place around 1955 between two of the upper-level managers of the Theoretical Division, one of whom said, ‘‘Look, what is this stuﬀ? This PIC method that theyÕre trying to publish a paper on will never work! There will never be computers big enough and fast enough to solve problems with all of the particles that you need in order to resolve whatÕs going on. I donÕt think we should publish this.’’ And the other person said, ‘‘Hey look, maybe itÕs not useful or wonÕt work all that well. But on the other hand, itÕs something that was tried, something that people will learn from. We probably should mark it down and itÕll go into the annals of things that were tried and didnÕt work very well but at least it will have been recorded.’’ IÕm happy that the second point of view prevailed and that the ideas found their way into print. Another type of opposition occurred in our interactions with editors of professional journals, and with scientists and engineers at various universities and industrial laboratories. One of the things we discovered in the 1950s and early 1960s was that there was a lot of suspicion about numerical techniques. Computers and the solutions you could calculate were said to be the playthings of rich laboratories. You couldnÕt learn very much unless you did studies analytically. The truth, of course, has turned out to lie in the complementary interaction between calculations and analysis, and in the validation of both through comparisons with experiments. I think some of these early suspicions were partially justiﬁed, although we found as time went by that many universities, industries, and governmental agencies began to solve problems of interest using numerical techniques. In the writing of reports and papers for publication we felt that it was important for readers who might wish to try our techniques that the writing be as complete as possible, even to the inclusion of the lengthy ﬁnite-diﬀerence equations. We were eager to have conﬁrmation from other investigators that the techniques could be used on the various computers that were rapidly becoming available to the world, and give meaningful results for solving actual problems that could arise in both engineering and scientiﬁc investigations. The part of the presentations that we found to be especially important was, and continues to be, the graphics. To absorb the meaning of vast amounts of computed data requires for all of us the pictorial depiction of results in ways that display the essence of what is going on. At ﬁrst our visual displays were obtained by the tried and true method of plotting by hand. Hundreds of PIC particles have been plotted by that tedious procedure. A tremendous improvement came with the availability of the Cal-Comp machines, which mechanically plotted the position of each particle with much greater speed and accuracy than we could accomplish by hand. But the greatest achievement of convenience and accuracy came with the capability of the computer itself to produce the desired graphics. Microﬁsch ﬁlms became our tool for this purpose, and many hours were spent with a Microﬁsch reader examining the graphs and pictures and trying to determine the best way to interpret the results of our calculations. The next improvement came with the capability for the computer to record the conﬁgurations of our ﬂow calculations on 35-mm ﬁlm. The result was our ﬁrst opportunity to make movies of the changing structure of the ﬂow! When I was told that COLOR depictions were soon to be possible, my ﬁrst reaction was to think that we already had all of the 420 F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 capability that we needed, and color would be an excessive luxury. How very wrong this old-fashioned viewpoint has turned out to be. One challenge that gave us some grief was how to develop good hidden-line routines. With all the current superb graphics packages available now itÕs amusing to remember how primitive our activities were in those days. 6. Early diﬃculties (1960, 1970) We had some exciting successes in the early days of the PIC, MAC, and ICE developments and applications, but not every numerical method that we tried worked out well. One example that we tried for compressible ﬂow was our experimentation with some ideas for eliminating the underlying mesh of computational cells, and putting all of the calculations onto the particles themselves. Each particle had its own mass, velocity, and internal energy, and the density associated with its location was to be obtained by examining the positions of the nearest neighbors [12,17,25]. We called the technique a Particle-and-Force (PAF) method because the dynamics of the system were determined on the basis of inter-particle forces derived in such a way as to represent the eﬀects of an appropriate equation of state for the material. A rather elaborate statistical theory was developed for the formulation of the forces. The idea is the forerunner of various techniques referred to, collectively, as Free Lagrangian. At our stage of the development, the PAF technique worked fairly well for a single material in two dimensions. Its main diﬃculties were seen in application to problems with two or more interacting materials. We also were hampered by the inability, at that time, for determining nearest-neighbor particles in a consistent way that was not excessive in computer resources. Bill Meixner, Bart Daly, and Eddie Welch joined with me in writing several reports on our ﬁndings, but we have not pursued this line of investigation any further, especially because of the brilliant accomplishments for this type of computations developed elsewhere. A second line of investigation resulted in what we call the Dynamics-of-Contours (DOC) method . The main feature of this approach is the automatic creation of a computational mesh of cells that is neither Lagrangian nor Eulerian, but is composed instead of cell boundaries that follow contours of constant density. This approach works ﬁne for one-dimensional problems but has signiﬁcant diﬃculty when applied to two-dimensional conﬁgurations subjected to strong distortions. Some of this diﬃculty was later removed by Norman Zabusky and other investigators, but the technique has not received much acceptance by potential users. 7. Multi-ﬁeld ﬂows (1971) Having worked with marker particles in computational cells it seemed like a natural extension for us to consider the transport of real, physical particles carried by a surrounding ﬂuid. The process is called multiﬁeld ﬂow. Multi-ﬁeld ﬂow refers, for example, to the motion of bubbles in oil, to raindrops in the air, to particles of sand in water, and to dust in the atmosphere. Much previous work had been done to develop the governing equations for this type of dynamics, based on straight-forward extensions of the Navier– Stokes equations for the dynamics of a purely homogeneous ﬂuid. The particles or bubbles and the ﬂuid are both represented by continuous ﬁelds that can interpenetrate. Local pressure equilibrium is assumed between the ﬁelds. As a result, the formulation contains one especially disturbing feature: it is mathematically ill-posed, which means that as you get to smaller and smaller scales of ﬂow, the solutions start wandering oﬀ in directions of their own, ultimately becoming arbitrarily far from each other, even for very small diﬀerences in the initial conditions. Another related diﬃculty arose in trying to represent the processes that take place on the scale of each individual particle, in particular the detailed exchanges of mass, momentum and energy. The concept of F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 421 continuous ﬁelds for both components of the materials means that these exchanges can only be described in a bulk constitutive fashion. As a result, there was no way that the dynamics could be expected to follow the formulation for scales smaller than the inter-particle spacing. Thus Tony Amsden and I built our solution technique for multi-ﬁeld ﬂow on the assumption that the only results that could be meaningful are those that occur at scales signiﬁcantly larger than the average value of that spacing, so that the problem of illposed behavior did not trouble us very much . We soon realized that any corrective procedure to remove the ill-posed nature of the equations would be ﬁne to incorporate as long as it did not interfere with the solutions at those scales for which the equations are physically meaningful. In particular, the existence of physical instability, which is related to the transition of the ﬂow to turbulent, needs to be preserved at all the scales that are relevant to the problem of interest. A little experimentation showed that the addition of a small amount of artiﬁcial viscosity to each ﬁeld works very well in curing the problem. To solve problems we used a straight-forward extension of the MAC-method ideas, employing the Newton–Raphson technique in each computational cycle to ﬁnd the pressure that brings the velocities to values that ensure overall conservation of mass. A more challenging problem was that of determining meaningful descriptions of the exchange processes. The drag of the interpenetrating ﬁelds upon each other was especially important to represent accurately. We postulated that the momentum exchange takes place by means of combined viscous drag and form drag, but were soon concerned that the formulation was based on the identiﬁcation of which was the dispersed ﬁeld, and which was the surrounding continuous ﬁeld . For widely separated spherical particles in water the identiﬁcation was not ambiguous, and the drag formula was reasonably good as long as there were no appreciable particle-to-particle interactions. For bubbles in water, however, the ﬂow circumstances could easily reverse the identiﬁcation to that of droplets in air. We developed a procedure that covers the behavior through the transition from one circumstance to the other, and it has worked reasonably well for many problems of interest to us. We have even applied the momentum-exchange equations to circumstances in which the entities (bubbles or droplets) are expected to deform appreciably during the dynamics, and we still got reasonably good results. The other exchange processes (of mass and energy) were especially important for problems in which there are phase transitions, as for example in steam-water ﬂows. Here was a crucial area for our developments, because of an assignment to consider the behavior of pressurized water reactors during a hypothetical blow-down in which you might have very hot pressurized water suddenly released to atmospheric pressure, forming bubbles and droplets and sending a potentially disastrous pressure pulse through the entire piping system of the reactor . Applications by Bart Daly and others in our Group to reactor safety analysis under the auspices of the US Nuclear Regulatory Agency were remarkably successful, and the initial studies proved so meaningful that a whole new Division of the Laboratory was created to extend our multi-ﬁeld technique to the calculation of the entire hot-water piping system of the reactor. Many other types of problems were also solved in our Group, including the transport of pollutants in the atmosphere around buildings in an urban setting , and the dynamics of fuel droplets in the cylinder of an internal combustion engine. These activities received substantial funding from agencies and industries outside the Laboratory. Many other examples are given in [64–67,81,99,100]. 8. Relativistic ﬂuid dynamics (1975) A completely diﬀerent type of application of our numerical techniques arose with problems of relativistic ﬂuid dynamics. Large accelerators at several Laboratories in the world were obtaining data on the eﬀects of atomic-nucleus collisions at very high energies for which the eﬀects of EinsteinÕs relativity are very important. To address the problem of interpreting these experimental results, Ray Nix from our nuclear 422 F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 physics Group and Fred Goldhaber from the State University of New York were discussing the possible applicability of a liquid-drop model. With their speciﬁcation of the appropriate relativistic ﬂuid dynamics equations, Tony Amsden and I set to work to modify the PIC method for applicability to such problems as the impact of a neon nucleus on the nucleus of a uranium atom. We put the uranium atom into the laboratory frame of the computational mesh, so that in the initial conditions the neon nucleus was appreciable shortened in the direction of its trajectory. One of the goals of the study was to obtain the distribution of fragments that propagated away from the point of collision, and the results of the calculations were remarkably successful in duplicating the experimental results. Six papers were published on our part of the project [68,69,71,72,74,76], and other investigators then took the ideas and have been continuing the studies with PIC-method calculations even to the present time. 9. The spread of numerical ﬂuid dynamics The development of numerical methods for solving material dynamics problems of a wide variety of kinds occupied my attention for many years. Up to about 1968, I felt conﬁdent in knowing about the dissemination and extensions of our techniques to most parts of the world, but at the same time began to realize the genius of other method developers in creating new techniques. Numerical ﬂuid dynamics had become an important and respectable science, and my directions of investigation began to change. 10. Turbulence (1967) We had long realized that turbulence occurs in many of the ﬂuid-ﬂow circumstances of interest to both scientiﬁc studies and engineering applications. Turbulence is like any of the ﬁne-scale processes that we have to represent in the behavior of a ﬂuid, whether it be a gas, a liquid, a plastic ﬂowing metal, or any other deformable material. It is like heat, which represents the eﬀects of molecular ﬂuctuations. Thermal energy describes part of the collective behavior of the molecules, but its transport does not require the determination of every molecular trajectory. The thermal energy per unit mass of the material is transported as though the molecular assemblage were continuous. Could turbulence behavior likewise be described by the evolution of some appropriately-chosen ﬁeld variables? The answer, of course, had been partially known for many years. Turbulence energy transport equations had been proposed by Rota in Germany and by P.Y. Chow in China, and had subsequently been generalized to transport the full Reynolds stress, but there was a big unresolved problem: Can we ﬁnd a transport equation for the determination of an appropriate length scale from which to calculate several required rate functions? These rates in the transport equation include the rates of return to isotropy and of decay. To describe the evolution of an appropriate length scale, Paul Nakayama and I introduced a transport equation for the rate of energy dissipation, epsilon [38,45]. Thus was born the K–epsilon model that has been used in our calculations and in computations performed by many other investigators throughout the world. A notable example is the work of B.E. Launder and D.B. Spalding. The model turned out to be very useful, and we tried solving many problems with it. The results often worked out nicely, but sometimes they were not as good as we wished, especially when the mean ﬂow of the ﬂuid was rapidly changing. Only somewhat later did we discover a spectral approach to the analysis that helped to extend the range of applicability. It was, however, a new and exciting sort of adventure and we enjoyed the challenge of testing the model by applying it to a wide variety of ﬂuid-ﬂow circumstances. One class of problems that we could not address with the ﬁrst version of these turbulence transport equations involved circumstances in which there are two or more ﬂuids with diﬀerent densities or a single ﬂuid with continuous variations of density. At the same time we were interested in the application of our F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 423 model to ﬂuids that are compressible, sometimes with strong shocks or rarefactions. To describe the eﬀects of variable density required the introduction of two new descriptive variables together with transport equations to describe their evolution. The two new descriptive quantities are the variance of the density ﬂuctuations and the mass ﬂux relative to a coordinate system in which the volume ﬂux is zero. The basic Navier–Stokes equations for the derivations are the same as we used for obtaining the constant-density K– epsilon model. By means of suitable manipulations, Rick Rauenzahn and Didier Besnard (visiting from France) joined with me in adding two more equations to describe the transport of the two new quantities [112,129]. Using the initials of the developers, the new technique came to be known as the BHR method. The ﬁrst application was to Rayleigh–Taylor instability, in which the mixing interpenetration of the materials across an interface is driven by a normal pressure gradient. Comparisons with abundant experimental data that had been obtained in British and Russian Laboratories showed good agreement for the growth of mixing-layer thickness even into the far non-linear regime. The results for the mean density proﬁle across the layer, however, were not satisfactory. The root of the problem lay in troubles associated with the calculation of the length scale, which was too small at the edges of the mixing layer. The missing element of the theory was demonstrated to lie in the rapid pressure-wave propagation to the edges of the layer. Using as guidance the linearized equations for low-amplitude Rayleigh–Taylor instability analysis, Michael Steinkamp and I saw that the missing element of the model was associated with the ‘‘non-local’’ spreading of the long-wave-length disturbances, which transport important acceleration to the ﬂuid beyond the edges of the mixing layer. With an integral spreading function for the Reynolds-stress source term the results were excellent. Many other unstable mixing processes have been investigated with considerable success. It should be noted that a turbulence transport model for variable-density ﬂuids, developed independently by scientists in the Soviet Union, is very similar to the BHR model. One of the interesting things that we realized at that time was that there are two major diﬀerences between turbulence and heat. With heat there is usually a rapid return to some kind of equilibrium after the ﬂuid is subjected to a disturbance. In addition, the length scale associated with molecular ﬂuctuations is essentially the same as the mean free path, for which the variations are easily calculated. With turbulence, the behavior is quite diﬀerent; there is a distribution of length scales that can vary signiﬁcantly with position and time, and the return to equilibrium or to any kind of self-similarity is very sluggish, taking place in times that are comparable to the overall elapsed time of interest. If you disturb a material with a rapid transient of drive, it will respond with the development of a strong departure from equilibrium, which persists for considerable time. To be more speciﬁc, we mean by non-equilibrium the distortion of the distribution of Reynolds-stress components from a ‘‘universal’’ self-similar form. The distribution of energy among all the wavelengths can look very strange, with humps and gaps in its structure. The idea of strong problem dependence on the distribution of energy across the size scales leads to the speculation that a spectral theory for the turbulence might be appropriate. It was felt that a close relation to the behavior of gas molecules could be found. For a gas with a mean free path that is comparable to the size of the macroscopic scales of the conﬁguration, the idea of a nearly universal distribution of molecular velocities (the Maxwell–Boltzmann distribution) was well known to be incorrect in circumstances with rapid changes in the dynamical drive of the system, or with large diﬀerences in the temperature and/or velocity of adjacent walls. For such circumstances it was necessary to transport a probability distribution function for the velocities. Could a similar analysis procedure be found for turbulence? We learned that exploration of this idea had taken place in France some years earlier, with results that seemed promising as a basis for developing the type of transport equations that are needed for this type of circumstances. To proceed with the derivation, Chuck Zemach joined with us in writing the Navier–Stokes equations for two diﬀerent points in space, so as to derive a way to examine the correlation of ﬂuctuations that could occur in the dynamics at those two positions . As the distance becomes large the correlations vanish. Thus a Fourier decomposition of the two-point equations with respect to the separation variable could be expected to describe the dynamics of the distribution function for such two-point correlations as those of 424 F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 velocity ﬂuctuations (the spectral form of the Reynolds stress). In many respects the result resembles the Boltzmann equation. In particular, appropriate moments of the equation over all of spectral space reduce to the K–epsilon equations if the assumption of a ‘‘universal’’ form for the spectral distribution is made. When that assumption is not valid, the spectral equations themselves can be solved to follow the evolution of the non-equilibrium behavior. To test the spectral equations Tim Clark and Patrick Spitz (visiting from France) extended them to include the presence of density variations, and we applied them to both pressure-driven and shear-driven interface instabilities. The buoyancy calculations, including the important non-locality of pressure-wave eﬀects in physical space, gave excellent results. The technique, described in several publications by Mike Steinkamp, Tim Clark and myself, has been called the SCH method [134,135]. For shear-driven instability, there was, however, a missing feature in the equations associated with a well-known process that continuously creates vortices of larger and larger sizes (called vortex pairing in a free shear layer). To represent this process, Peter Wilson and I found that a modiﬁed version of the Kraichnan-Spiegel non-locality formula in spectral space could be made to work very well. With our spectral–transport model, Peter Wilson and I decided to test the equations with one of the classic non-equilibrium (non-self-similar) problems, namely that of the combined pressure-driven AND shear-driven instabilities at an interface between two ﬂuids. The calculations for each separately had given excellent results for the development of the self-similar spectral structures. With our new modiﬁcations, the calculations showed that the presence of both shear and buoyancy produces a continuing non-equilibrium in which the spectral structure has ever-changing features that continue until the occurrence of ultimate late-time dominance of the buoyancy eﬀects. Perhaps the most useful of the results of the spectral derivations is that they furnish a procedure for solving problems numerically for which the K–epsilon model and its variants are not capable of capturing important eﬀects of spectral non-equilibrium. They have also furnished for the ﬁrst time a procedure by which the evolution equation for epsilon could be derived, instead of simply postulated on the basis of dimensional and geometric arguments. 11. Some applications of numerical methods In the development of new solution techniques for ﬂuid dynamics we needed to perform many types of test problems in order to validate the relevance and accuracy of the procedures. A few of these test problems are described above. At the same time, however, we undertook a variety of large-scale applications, ﬁrst because they were interesting and relevant ways to increase our conﬁdence in these new techniques, and second for the pragmatic purpose of developing sources of funding for our work. The Laboratory, along with many other scientiﬁc institutions in this country, found that the federal support for research was dwindling in the late 1960s, so that we needed to make contacts with various other government agencies that were still able to support projects of interest to their speciﬁc needs. We had computational tools that were new and powerful, and soon found that there was great interest in other agencies for us to apply our capabilities to problems of relevance to their requirements. A small number of examples are mentioned here to indicate the scope of our activities. The reactor-safety analysis program is one of these projects, which is described above. Another, undertaken under the auspices of the US Environmental Protection Agency, was the analysis of pollution transport by winds through urban areas with streets and tall buildings. Bob Hotchkiss, Tony Hirt, Larry Cook and I used an extension of the MAC method with the inclusion of heat transport for representing buoyancy eﬀects, and a particulate transport equation to describe the combined eﬀects of gravity-driven settling, drag with the air, and turbulent diﬀusion. Sources of the pollution included the emission from automobiles and the eﬄuence of material from rooftop vents. The calculations were applied to various F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 425 conﬁgurations of buildings and circumstances of wind direction and speed . More particularly, we examined the pollution transport in the Broadway Street Canyon in downtown St. Louis, Missouri. The calculations worked so well that many research institutions built on our methodology to include complex terrain and many other factors contributing to the reality of speciﬁc pollution-transport circumstances. Geothermal energy extraction was another challenge that we addressed, using our multi-ﬁeld solution techniques. The particular circumstances of interest to one of the experimental Divisions of our Laboratory involved the hot dry rocks that lie buried in the heart of the mountains just to the west of Los Alamos. The technique employed in the ﬁeld requires two deep holes to be drilled into the rocks. Cold water is pumped down into one of the holes, is heated by the rocks, and extracted through the other hole. Success requires a connection between the holes, and especially that the cold water fracture the rocks in order that the ﬂow region would continuously enlarge for many years of energy extraction. Bill Pracht and I used a multi-ﬁeld approach to study the problem, and the results we obtained continued to inﬂuence the program for many years . Another thing that was of considerable interest in the early 1970s was the behavior of tornadoes. At that time there had been a major tornado in the city of Lubbock, Texas. This was an incredibly devastating storm, which did not behave like an ordinary tornado; there was something strange and diﬀerent about it. To investigate what had happened, the US National Science Foundation gave us a grant to apply our numerical techniques to the analysis. Lee Stein and I calculated the dynamics of converging ﬂow with various magnitudes of preexisting vorticity, and found that we could distinguish between a single-cell tornado and a double-cell tornado. In the ﬁrst type the potential for damage is strong; while in the second type the potential for damage is enormously greater. Within the single-cell structure the central updraft extends upward into the spawning clouds, whereas within the double-cell tornado the centrifugal acceleration is so great that the air is thrown outward, and above a vertical stagnation point the velocity reverses and air travels downward from the clouds. The violence of the rotation can lead to the formation of a set of smaller tornadoes traveling around the central one, whose paths had been observed in damage patterns after the storm was over, and also has been observed in dry tornadoes, often called dust devils. For the milder type of tornado we were able to compare our results with those that had been obtained in a laboratory vortex chamber at another institution, with results that lent signiﬁcant conﬁdence to the validity of our calculations [61,63]. A very diﬀerent type of study that we did was to investigate the dynamics of a biological cell, in which the internal circulation can give the cell the ability to crawl along a surface or to send out ﬁngers that are called pseudopods. The work required the collaboration of Micah Dembo, a biologist in our Laboratory, because many of the relevant aspects of the internal structure of the cell were completely unknown to Mat Maltrud and me. We were very excited about the calculations, and enjoyed the publicity that came from publication of the results in several journal articles [106,108]. A fascinating project that I especially enjoyed was the development of a large computer code to describe the behavior of a wildﬁre propagating through a forest. The conﬁguration that Rod Linn and I wanted to represent is fully three-dimensional; the forest is situated on rugged terrain consisting of deep canyons, riverbeds, ﬂat-topped mesas, and various other relevant geological formations that might be present. The challenge is to incorporate into the code complex physical processes over a large range of scales, along with the chemistry of hydrocarbon release and combustion, and with boundary conditions representing realistic weather conditions such as the time-varying wind speed, intensity and scale of turbulence, and humidity. Self-determination of propagation rates, transition to crowning, automatic choice of the most likely burn path, and the eﬀects of spotting are among the answers that we wanted to extract from the calculations. It was clear from the start of this project that no calculation is ever going to be able to predict the exact outcome of a particular wildﬁre; the best that can be expected is a meaningful description of the relative probabilities for the diﬀerent possible scenarios to develop. The practical value of these calculations is to complement the excellent rapid-scoping techniques that have been developed by the US Forest Service, and 426 F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 by other ﬁre-mitigation agencies for estimating the behavior of wildﬁres. Some of the goals include the development of urban-planning strategies to minimize potential threats to new housing developments, the training of ﬁre ﬁghters, and the planning of controlled burns in forest regions that are considered to be especially susceptible to disastrous wildﬁre activity. Combining our developments with the atmospheric code of Jon Reisner, we tested the calculations by comparison with well-documented actual ﬁres in California, Colorado, New Mexico, and Florida, with excellent results [136,141]. 12. Scale bridging In recent years, I have become strongly interested in the potential for analytical–computational studies of systems with many degrees of freedom (like our turbulence studies using Reynolds-decomposition techniques), in which a probability distribution function (PDF) can be evolved by means of equations based on the Liouville principles of conservation in a generalized phase space. We have applied these ideas to a variety of problems. One of them, with Mike Steinzig, is to the development of grains in the solidiﬁcation of a binary alloy, for which the mean grain size and distribution of sizes need to be calculated as functions of cooling rate and relative proportions of the two metals in the alloy . Another, with Eric Harstad, is to the behavior of a polymer molecule during the non-equilibrium dynamics taking place at high strain rates, for which the Liouville approach resembles that of the techniques used to study Brownian motion . The response of a polymeric foam to mechanical deformations has also been studied by Mark Schraad and me . We follow the evolution of a probability distribution function for the variations of strain and stress in a random distribution of cell sizes. To describe the eﬀects of strain-rate variations in the foam we include a version of the two-ﬁeld equations to represent the movement of air through the matrix of cells. 13. Social problems An interesting question arose from our investigations of scale-bridging techniques using PDF transport: If we can represent the combined behavior of individual molecules, and if we can represent the bulk eﬀects of turbulence eddies, can we go one step further and model the collective behavior of groups of people? To answer this question we worked on the development of model equations that we thought could be relevant. I had two summer students who worked with me, Don Sandoval and Kent Genin, and we developed essentially the same technique that works well for molecular dynamics and for turbulence dynamics, but involves one important new factor, and that is the presence of what we call soft variables . By soft variables, we mean properties in addition to such deterministic quantities as the position and the velocity of the individual person, namely the thoughts associated with excitement and fear. A major challenge, of course, is to quantify the intensity of these emotional responses and to relate them to the behavior of the individual. In a rampaging mob, for example, the descriptions can be expressed in terms of either running towards something if the personÕs excitement is great enough or running away from something if his fear is great enough. We developed a relatively simple model and published a report about our application to the behavior of two platoons of soldiers in battle with one another under a variety of diﬀerent circumstances. Now, almost 15 years later, I am ﬁnding that these types of models are starting to be of increasing interest to people, both in questions of international politics and economics, as well as the behavior of mobs of people, and most recently in issues of Homeland Defense. I predict that social dynamics (both slow and fast) with soft variables will become a major ﬁeld of computational activity in the coming decades. F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 427 14. Additional interests When I ﬁrst came here it seemed important to do things besides science. My family was and continues to be a major source of delight. To my wife, Patricia, I am especially grateful for more than 50 years of loving friendship and support. In my spare time, paleontology became a signiﬁcant activity for quite a while and Professor Patrick Sutherland from the University of Oklahoma and I wrote a major treatise on the fossil brachiopods of the upper Carboniferous of north-central New Mexico. I have also developed a lot of interest in the Indian pottery of the southwestern part of the United States. The Native Americans of this area are famous for manufacturing pottery using the same techniques developed many centuries ago. I have written six books on the subject and am working on three more, two of them already at the publisher; these are described in Section II of the Bibliography [152–174]. Also, since 1968, it has been both a challenge and a continuing source of recreation to paint well over a thousand pictures of the Indian pottery and of many other scenes that are typical of the activities of southwestern Indians, cowboys, and motorcycle riders. Acknowledgements I have lived a lucky life, being in the right place at the time when something new and exciting was happening. The activities have succeeded only through a combination of pure luck and pure joy, and especially as a result of having wonderful people to work with. Indeed each of the projects described in this paper could only have succeeded with the collaboration of many remarkable associates. The people with whom I published reports and papers during the pioneering years for each project are listed in the Bibliography. To each of my co-authors I shall always be more grateful than words can ever express. There have been many other people in our Group who have made marvelous progress in the ﬁeld of numerical ﬂuid dynamics and also of continuum dynamics in general, with whom I have interacted through many inspiring technical discussions. I am especially grateful to Frank Addessio, Tony Amsden, Didier Besnard, Jerry Brackbill, Dan Butler, Tim Clark, Tom Cook, Bart Daly, Ruth Demuth, Martha Evans, Jake Fromm, Dick Gentry, Rob Gore, Eric Harstad, Tony Hirt, John Hopson, Marty Horn, Bob Hotchkiss, Norm Johnson, Bucky Kashiwa, Dana Knoll, Rod Linn, Mat Maltrud, Len Margolin, Ron Moses, Paul Nakayama, Peter OÕRourke, Bill Pracht, Rick Rauenzahn, Murray Rudman, Hans Ruppel, Manjit Sahota, Don Sandoval, Mark Schraad, Mike Steinkamp, Mike Steinzig, Brian VanderHeyden, Eddie Welch, Peter Wilson, Chuck Zemach and Duan Zhang for many years of friendship, support, and numerous stimulating discussions. The publications listed below show only a part of the work accomplished by our Group. Virtually every new development saw many years of applications to problems, accomplished by later Group members and by researchers at many other places around the world. To all of them I wish to express my gratitude for their interest and for the wonderful feedback that they have expressed. In addition, there have been literally hundreds of university students who have taken advantage of this LaboratoryÕs student programs, and have come to Los Alamos for summer activities in our Group. Fifteen of them have stayed here for two or three years to work with me on the research and writing of their Ph.D. dissertations for nearly as many universities. All of them have brought their youthful enthusiasm to add immensely to the excitement of our work, and to the accomplishment of our missions. References Extensive listings of publications resulting from some of our early activities are contained in the papers: Francis H. Harlow, ‘‘PIC and its Progeny’’, Computer Phys. Com., Volume 48, pages 1–10 (1988). Margaret Findley, ‘‘Bibliography of Group T-3, 1955-1991’’, Los Alamos National Laboratory Report Number LA-12377-MS, September, 1992. 428 F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 My personal list of publications includes the following: I. Physics and Mathematics                                   F.H. Harlow, X-ray measurements of the cyclotron rf dee voltage, Rev. Sci. Instr. 23 (September) (1952) 501. F.H. Harlow, B.A. Jacobson, Nucleon isobars in intermediate coupling, Phys. Rev. 93 (January) (1954) 333. F.H. Harlow, Hydrodynamics problem involving large ﬂuid distortions, J. Assoc. Compos. Mach. 4 (April) (1957) 137. M.W. Evans, F.H. Harlow, The particle-in-cell method for hydrodynamic calculations, Los Alamos Scientiﬁc Laboratory report LA-2139 (1957). See also A.A. Amsden, The particle-in-cell method for the calculation of the dynamics of compressible ﬂuids, Los Alamos Scientiﬁc Laboratory report LA-3466, 1966. M.W. Evans, F.H. Harlow, Calculation of supersonic ﬂow past an axially symmetric cylinder, J. Aerosol Sci. 25 (April) (1958) 269. F.H. Harlow, B.D. Meixner, The diﬀusion of radiation, Los Alamos Scientiﬁc Laboratory report LA-2196, 1958. M.W. Evans, F.H. Harlow, Calculation of unsteady supersonic ﬂow past a circular cylinder, A.R.S. J. 29 (January) (1959) 46. F.H. Harlow, D.O. Dickman, Numerical study of the motions of variously-shaped slabs accelerated by a hot gas, Los Alamos Scientiﬁc Laboratory report LA-2256. F.H. Harlow, D.O. Dickman, D.E. Harris Jr., R.E. Martin, two-dimensional hydrodynamics calculations, Los Alamos Scientiﬁc Laboratory report LA-2301, 1959. F.H. Harlow, Stability of diﬀerence equations, Los Alamos Scientiﬁc Laboratory report LA-2452-MS, 1960. F.H. Harlow, Dynamics of compressible ﬂuids, Los Alamos Scientiﬁc Laboratory report LA-2412, 1960. F.H. Harlow, B.D. Meixner, The particle-and-force computing method for ﬂuid dynamics, Los Alamos Scientiﬁc Laboratory report LA-2567-MS, 1961. F.H. Harlow, B.D. Meixner, Motion of a viscous compressible gas adjacent to a sliding plate, Phys. Fluids 4 (October) (1961) 1202. F.H. Harlow, R.E. Duﬀ, C.W. Hirt, Eﬀects of diﬀusion on interface instability between gases, Phys. Fluids 5 (April) (1962) 417. F.H. Harlow, M.W. Evans, B.D. Meixner, Interaction of shock or rarefaction with a bubble, Phys. Fluids 5 (June) (1962) 651. F.H. Harlow, The particle-in-cell method for numerical solutions of problems in ﬂuid dynamics, in: Invited Lecture at Symp. on Experimental Arithmetic, American Mathematical Society, Chicago, April 12–14, 1962. Also published in Proc. Symp. Appl. Math., vol. 15, 1963, p. 269. F.H. Harlow, Theory of correspondence between ﬂuid dynamics and particle-and-force models, Los Alamos Scientiﬁc Laboratory report LA-2806, 1963. F.H. Harlow, J.E. Fromm, Numerical solution of the problem of vortex street development, Phys. Fluids 6 (July 1963) 975. Republished in the AlAA Series of Selected Reprints, vol. 4, in: C.K. Chu (Ed.), Computational Fluid Dynamics, Reprinted in Selected Papers in Physics, vol. VI, The Physical Society of Japan, Tokyo, 1971. F.H. Harlow, The particle-in-cell computing method for ﬂuid dynamics, Methods Comput. Phys. 3 (1964) 319–343. A.A. Amsden, F.H. Harlow, Slip instability, Phys. Fluids 7 (March) (1964) 327–334. R.A. Gentry, F.H. Harlow, R.E. Martin, Computer experiments for molecular dynamics problems, Methods Comput. Phys. 4 (1964) 211–245. F.H. Harlow, J.E. Fromm, Dynamics and heat transfer in the von karman wake of a rectangular cylinder, Phys. Fluids 7 (August) (1964) 1147–1156. F.H. Harlow, J.E. Fromm, Computer experiments in ﬂuid dynamics, Sci. Am. 212 (3) (1965) 104–110, Reprinted in Scientiﬁc American Reader on Computers and Computation, Joseph Weizenbaum and Robert Fenichel of MIT, Ed. (1971). F.H. Harlow, Numerical ﬂuid dynamics, Am. Math. Mon. 72 (2) (1965) 84–91. B.J. Daly, F.H. Harlow, J.E. Welch, E.E. Sanmann, E. Wilson, Numerical ﬂuid dynamics using the particle-and-force method, Los Alamos Scientiﬁc Laboratory report LA-3144. A.A. Amsden, F.H. Harlow, Numerical calculation of supersonic wake ﬂow, AIAA J. 3 (November) (1965) 2081–2086, Republished in the AlAA Series of Selected Reprints, vol. 4, Computational Fluid Dynamics, C.K. Chu, Ed. F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible ﬂow of ﬂuid with free surface, Phys. Fluids 8 (December) (1965) 2182–2189, Reprinted in Selected Papers in Physics, vol. VI, The Physical Society of Japan, Tokyo, 1971. F.H. Harlow, J.P. Shannon, J.E. Welch, Liquid waves by computer, Science 149 (September) (1965) 1092–1093. F.H. Harlow, J.E. Welch, Numerical study of large amplitude free surface motions, Phys. Fluids 9 (May) (1966) 842–851. F.H. Harlow, J.E. Welch, J.P. Shannon, B.J. Daly, The MAC method, Los Alamos Scientiﬁc Laboratory report LA-3425, March 1966. F.H. Harlow, P.I. Nakayama, Simulating ﬂuid turbulence, Sci. J. 3 (9) (1967) 74. F.H. Harlow, W.E. Pracht, Formation and penetration of high speed collapse jets, Phys. Fluids 9 (October) (1966) 1951–1959. F.H. Harlow, J.P. Shannon, J.E. Welch, Un Calculateur Qui Fait Des Vagues, Sciences 7 (41) (1966) 14–17. F.H. Harlow, P.J. Nakayama, Turbulence transport equations, Phys. Fluids 10 (November) (1967) 2323. F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 429  F.H. Harlow, J.P. Shannon, The splash of a liquid drop, J. Appl. Phys. 38 (September) (1967) 3855.  F.H. Harlow, C.W. Hirt, A general corrective procedure for the numerical solution of initial-value problems, J. Comput. Phys. 2 (1968) 114.  F.H. Harlow, J.P. Shannon, Distortion of a splashing liquid drop, Science 157 (August) (1967) 547–550.  F.H. Harlow, P.I. Nakayama, Transport of turbulence energy decay rate, Los Alamos Scientiﬁc Laboratory report LA-3854, 1968.  A.A. Amsden, F.H. Harlow, Transport of turbulence in numerical ﬂuid dynamics, J. Comput. Phys. 3 (1968) 94.  F.H. Harlow, A.A. Amsden, Numerical calculation of almost incompressible ﬂow, J. Comput. Phys. 3 (1968) 80.  F.H. Harlow, Transport of anisotropic or low intensity turbulence, Los Alamos Scientiﬁc Laboratory report LA-3947, 1968.  F.H. Harlow, C.W. Hirt, Generalized transport theory of anisotropic turbulence, Los Alamos Scientiﬁc Laboratory report LA4086, 1969.  F.H. Harlow, Computer solutions in continuum mechanics, in: S. Fembach, A. Taub (Eds.), Computers and Their Role in the Physical Sciences, Gordon and Breach Science, Publishers, New York, 1970.  F.H. Harlow, A.A. Amsden, Fluid dynamics – an introductory text, Los Alamos Scientiﬁc Laboratory report LA-4100, 1969.  B.J. Daly, F.H. Harlow, Transport equations in turbulence, Phys. Fluids 13 (1970) 2634.  F.H. Harlow, N.C. Romero, Turbulence distortion in a nonuniform tunnel, Los Alamos Scientiﬁc Laboratory report LA-4247, 1969.  F.H. Harlow, Numerical methods for ﬂuid dynamics – an annotated bibliography, Los Alamos Scientiﬁc Laboratory report LA4281, 1969.  F.H. Harlow, A.A. Amsden, A numerical ﬂuid dynamics calculation method for all ﬂow speeds, J. Comput. Phys. 8 (1971) 197– 213.  F.H. Harlow, A.A. Amsden, A simpliﬁed MAC technique for incompressible ﬂuid ﬂow calculations, J. Comput. Phys. 6 (1970) 322–325.  A.A. Amsden, F.H. Harlow, The SMAC method – a numerical technique for calculating incompressible ﬂuid ﬂows, Los Alamos Scientiﬁc Laboratory report LA-4370, 1970.  F.H. Harlow, Contour dynamics for numerical ﬂuid ﬂow calculations, J: Comput. Phys. 8 (1971) 214–229.  B.J. Daly, F.H. Harlow, Inclusion of turbulence eﬀects in numerical ﬂuid dynamics, in: Proc. of the Second Intern. Conf. on Numerical Methods in Fluid Dynamics, Univ. of Calif., Berkeley, CA, Sept. 15–19, 1970.  F.H. Harlow, A.A. Amsden, C.W. Hirt, Numerical calculation of ﬂuid ﬂows at arbitrary Mach number, in: Proc. of the Second Intern. Conf. on Numerical Methods in Fluid Dynamics, Univ. of Calif., Berkeley, CA, Sept. 15–19, 1970.  F.H. Harlow, A.A. Amsden, Fluid Dynamics, second ed., Los Alamos Scientiﬁc Laboratory Monograph, LA-4700, 1971.  F.H. Harlow, C.W. Hirt, Recent extensions to Eulerian methods for numerical ﬂuid dynamics, J. Comput. Math. Math. Phys. USSR 12 (3) (1972) 656–672.  F.H. Harlow, W.E. Pracht, A theoretical study of geothermal energy extraction, J. Geophys. Res. 77 (1972) 7038.  F.H. Harlow (Ed.), Turbulence Transport Modeling, AIAA Selected Reprints Series, vol. XIV, New York, 1973.  F.H. Harlow (Ed.), Computer Fluid Dynamics – Recent Advances, AIAA Selected Reprints Series, vol. XV, New York, 1973.  B.J. Daly, F.H. Harlow, Reply to Comments by S. Corrsin, Phys. Fluids 16 (1973) 158.  R.S. Hotchkiss, F.H. Harlow, Air pollution transport in street canyons, EPA-R473-029, June 1973.  F.H. Harlow, L.R. Stein, Structural analysis of tornado-like vortices, J. Atmos. Sci. 31 (1974) 2081.  F.H. Harlow, A.A. Amsden, Multi-ﬂuid ﬂow calculations at all Mach numbers, J. Comput. Phys. 16 (1974) 1.  L.R. Stein, F.H. Harlow, Numerical solution of the ﬂow structure in tornado-like vortices, Los Alamos Scientiﬁc Laboratory report LA-57l3-MS, 1974.  F.H. Harlow, A.A. Amsden, Numerical calculation of multiphase ﬂuid ﬂow, J. Comput. Phys. 17 (1975) 19.  A.A. Amsden, F.H. Harlow, KACHINA: an Eulerian Computer Program for multiﬁeld ﬂuid ﬂows, Los Alamos Scientiﬁc Laboratory report LA-5680, 1975.  F.H. Harlow, A.A. Amsden, Flow of interpenetrating material phases, J. Comput. Phys. 18 (1975) 440.  J.R. Travis, F.H. Harlow, A.A. Amsden, Numerical calculation of two-phase ﬂows, Los Alamos Scientiﬁc Laboratory report LA-5942-MS (1975); also published in Nucl. Sci. Eng. 61 (1976) 1.  F.H. Harlow, A.A. Amsden, J.R. Nix, Relativistic ﬂuid dynamics calculations with the particle-in-cell technique, J. Comput. Phys. 20 (1976) 119.  A.A. Amsden, G.F. Bertsch, F.H. Harlow, J.R. Nix, Relativistic hydrodynamic theory of heavy-ion collisions, Phys. Rev. Lett. 35 (1975) 905.  J.W. Taylor, F.H. Harlow, A.A. Amsden, Dynamical plastic instabilities in stretching plates and shells, J. Appl. Mech. 45 (1978) 105.  A.A. Amsden, F.H. Harlow, J.R. Nix, Relativistic nuclear ﬂuid dynamics, Phys. Rev. C 15 (1977) 2059, Reprinted in: Kozi Nakai, Isao Tanihata, Koichi Tazaki (Eds.), High-Energy Nucleus–Nucleus Collisions, Selected Papers in Physics, The Physics Society of Japan, 1982. 430 F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433  A.A. Amsden, A.S. Goldhaber, F.H. Harlow, J.R. Nix, Relativistic two-ﬂuid model of nucleus-nucleus collisions, Phys. Rev. C 17 (1978) 2080.  F.H. Harlow, T.L. Cook, three-dimensional dynamics of rotating protostars, in: 149th meeting of the American Astronomical Society, Honolulu, January 1977.  A.A. Amsden, J.N. Ginocchio, F.H. Harlow, J.R. Nix, M. Danose, E.C. Hlabert, R.K. Smith Jr., Comparison of macroscopic and microscopic calculations of high-energy 20Ne +238U collisions, Phys. Rev. Lett. 38 (1977) 1055.  A.A. Amsden, T.D. Butler, F.H. Harlow, Numerical study of downcomer ﬂow dynamics, Los Alamos Scientiﬁc Laboratory report LA-NUREG-6797-SR, 1977.  A.A. Amsden, A.S. Goldhaber, F.H. Harlow, P. Moller, J.R. Nix, A.J. Sierk, Hydrodynamics calculations of heavy-ion collisions, in: Proc. of the Symp. on Macroscopic Features of Heavy-Ion Collisions and Pre-Equilibrium Processes, Hakone, Japan, September 2–3, 1977.  A.A. Amsden, F.H. Harlow, K-TIF: a two-ﬁeld computer program for downcomer ﬂow dynamics, Los Alamos Scientiﬁc Laboratory report LA-NUREG-6994, 1977.  T.L. Cook, F.H. Harlow, Three-dimensional dynamics of protostellar evolution, Astrophys. J. 225 (1978) 1005.  T.L. Cook, R.B. Demuth, F.H. Harlow, Numerical calculation of high-speed multiphase ﬂow, in: 31st Meeting of Fluid Dynamics Division, A.P.S., 19–21 November 1978.  T.L. Cook, R.B. Demuth, F.H. Harlow, Multiphase interpenetration of shocked materials, Los Alamos Scientiﬁc Laboratory report LA-7578, 1979.  T.L. Cook, R.B. Demuth, F.H. Harlow, PIC calculations of multiphase ﬂow, J. Comput. Phys. 41 (1981) 51.  J. Daly, F.H. Harlow, Scaling and constitutive relationships in downcomer modeling, Los Alamos Scientiﬁc Laboratory report LA-7610, 1979.  R.B. Demuth, F.H. Harlow, Thermal fracture eﬀects in geothermal energy extraction, Los Alamos Scientiﬁc Laboratory report LA-7963, 1979.  R.B. Demuth, F.H. Harlow, Geothermal energy enhancement by thermal fracture, Los Alamos Scientiﬁc Laboratory report LA8428, 1980.  F.H. Harlow, Mental dynamics, Los Alamos Scientiﬁc Laboratory report LA-7946-MS, 1979.  J. Daly, F.H. Harlow, Numerical study of condensation in concurrent stratiﬁed ﬂows, Los Alamos Scientiﬁc Laboratory report LA-8077, NUREG/CR-II08, 1979.  F.H. Harlow, B.A. Kashiwa, An investigation of simultaneous heat and mass transfer in subbituminous coal, in: Proc. of the 15th Intersociety Energy Conversion Engineering Conf., Seattle, WA, August 1980.  F.H. Harlow, B.A. Kashiwa, Hot gas drying calculations for a coal seam channel, in: Proc. of the Sixth Underground Coal Conversion Symp., Afton, OK, July 13–17, 1980.  J. Daly, F.H. Harlow, A model of countercurrent steam-water ﬂow in large horizontal pipes, Nucl. Sci. Eng. 77 (1981) 273.  A.A. Amsden, J.K. Dienes, F.H. Harlow, H.M. Ruppel, three-dimensional analysis of ﬂuid-structure interactions during blowdown of an Upper Plenum, Los Alamos Scientiﬁc Laboratory report LA-8754-MS, 1981.  F.H. Harlow, H.M. Ruppel, Propagation of a liquid–liquid explosion, Los Alamos National Laboratory report LA-897l-MS, 1981.  B.A. Kashiwa, F.H. Harlow, R.B. Demuth, T.J. Ruppel, Steam–water jet analysis, EPRI Research Project RP-1324-4, Final Report, May 1984.  T.L. Cook, F.H. Harlow, J.R. Travis, T.J. Bartel, C.E. Tyner, Heat and mass transfer in porous media, in: Proc. of the Fourteenth Oil Shale Symposium, Colorado School of Mines, Golden, Colorado, 22–24 April 1981.  J.R. Travis, F.H. Harlow, M.D. Torrey, Modeling of H2 migration in LWR containments, in: Proc. of the Ninth Water Reactor Safety Information Meeting, October 26–30, 1981, Gaithersburg, MD.  M. Dembo, F.H. Harlow, W. Alt, The biophysics of cell surface motility, in: A. Perelson, C. De Lisi, F. Weigel (Eds.), Cell Surface Dynamics: Concepts and Models, Marcel Decker, New York, 1984.  B.J. Daly, F.H. Harlow, Droplet de-entrainment and re-entrainment in the SCTF Upper Plenum – A Status Report, 2D/3D Program Technical Note, N.R.C., 1982 (LA-UR- 82-836).  F.H. Harlow, Numerical ﬂuid dynamics, Videotape Course of 15 lectures, Los Alamos National Laboratory, January–February 1982.  F.H. Harlow, N. Metropolis, Computing and computers at Los Alamos, Los Alamos Sci. (7) (1983) 132–141.  T.L. Cook, F.H. Harlow, Virtual mass in multiphase ﬂow, Int. J. Multiphase Flow 10 (1984) 691.  T.L. Cook, F.H. Harlow, Vortices in bubbly two-phase ﬂow, Int. J. Multiphase Flow 12 (1986) 35.  T.L. Cook, F.H. Harlow, VORT: a computer code for bubbly two-phase ﬂow, LO S 1. Alamos National Laboratory report LAIO021-MS, 1984.  F.H. Harlow, Super problems for supercomputers, in: N. Metropolis, D.H. Sharp, W.J. Worlton, K.R. Ames (Eds.), Frontiers of Super-Computing, University of California Press, Berkeley, 1986.  D. Besnard, F.H. Harlow, Turbulence in two-ﬁeld incompressible ﬂow, Los Alamos National Laboratory report LA-10187-MS, 1985. F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 431  D. Besnard, F.H. Harlow, Non-spherical particles in two-phase ﬂow, Int. J. Multiphase Flow 12 (1986) 891.  M. Raju, J. Phillips, F.H. Harlow (Eds.), Creativity in Science, Proc. of a Los Alamos Symp., 13–14 August 1984; Los Alamos National Laboratory report, LA-10490-C, 1985.  M. Dembo, F.H. Harlow, Cell motion, contractile networks, and the physics of interpenetrating reactive ﬂow, Biophys. J. 50 (1986) 109.  F.H. Harlow, D. Besnard, Well-posed two-phase ﬂow equations with turbulence transport, Lett. Math. Phys. 10 (1985) 161.  M. Dembo, M. Maltrud, F.H. Harlow, Numerical studies of unreactive contractile networks, Biophys. J. 50 (1986) 123.  F.H. Harlow, D.L. Sandoval, Mathematical model for mental activity, in: Mathematical Modeling of Biological Ensembles, Los Alamos National Laboratory report LA-10765-MS, 1986.  F.H. Harlow, D.L. Sandoval, Human collective dynamics: the mathematical modeling of mobs, in: Mathematical Modeling of Biological Ensembles, Los Alamos National Laboratory report LA-10765-MS, 1986.  F.H. Harlow, D.L. Sandoval, H.M. Ruppel, Mathematical modeling of biological evolution: a framework for the investigation of unanswered questions, in: Mathematical Modeling of Biological Ensembles, Los Alamos National Laboratory report LA10765-MS, 1986.  D. Besnard, F.H. Harlow, R.M. Rauenzahn, Conservation and transport properties of turbulence with large density variations, Los Alamos National Laboratory report LA-10911-MS, 1987.  D. Besnard, F.H. Harlow, Sources of turbulence in ﬂuid ﬂow, Los Alamos National Laboratory, Institutional Supporting Research and Development, Annual Report LA-10600, 1985.  D. Besnard, F.H. Harlow, N.L. Johnson, R. Rauenzahn, J. Wolfe, Turbulence transport, Los Alamos Sci. (Ulam Memorial Issue), 1987.  D. Besnard, J.F. Haas, M. Bonnet, A. Froger, S. Gauthier, B. Silt, F.H. Harlow, Comparison of two models of Rayleigh–Taylor induced turbulent mixing, in: Proc. of the Los Alamos/Limeil Conference on Numerical Methods in High-Temperature Physics, Los Alamos National Laboratory report LA-l1342-C, 1988.  D. Besnard, R. Rauenzahn, F.H. Harlow, Turbulence theory for material mixtures, in: Proc. of the Los Alamos/Limeil Conference on Numerical Methods in High-Temperature Physics, Los Alamos National Laboratory report LA-11342-C, 1988.  D.C. Besnard, R.M. Rauenzahn, F.H. Harlow, On the modeling of turbulence for material mixtures, in: G. de Vahl Davis, C. Fletcher (Eds.), Computational Fluid Dynamics, North-Holland, Amsterdam, 1988.  D. Besnard, F.H. Harlow, Un Modele de Turbulence dans les Melanges, part 2: Transport de la Turbulence et Establissement des Melanges, Commissariat a lÕEnergie Atomique, France, Special Report, 1987.  F.H. Harlow, Turbulence, in: Proceedings of Defense Nuclear Agency Numerical Methods Symposium, SRI International, April 28–29, 1987.  F.H. Harlow, PIC and its progeny, in: Proceedings of Los Alamos Workshop on Particle Methods in Fluid Dynamics and Plasma Physics, Los Alamos National Laboratory, April 13–15, 1987; Comput. Phys. Comm. 48 (1) (1988); also reprinted in book form, Particle Methods in Fluid Dynamics and Plasma Physics, North-Holland, Amsterdam, 1988.  D.L. Sandoval, F.H. Harlow, K.E. Genin, Human collective dynamics: two groups in Adversarial Encounter, Los Alamos National Laboratory report LA-11247-MS, 1988.  D. Besnard, F.H. Harlow, Turbulence in multiphase ﬂow, Int. J. Multiphase Flow 14 (1988) 679.  D. Besnard, F. Harlow, R. Rauenzahn, C. Zemach, Spectral model for inhomogeneous turbulence, in: XIXth Biennial Symp. on Fluid Mechanics, Polish Academy of Sciences, Kozubnik, Poland, September 3–8, 1989.  D. Besnard, F. Harlow, R. Rauenzahn, C. Zemach, Spectral transport model for turbulence, Los Alamos National Laboratory report LA-11821-MS (1990); also published in Theoret. Comput. Fluid Dyn. 8 (1996) 1–35.  D. Besnard, F. Harlow, R. Rauenzahn, C. Zemach, Un Modele Spectral de Transport de la Turbulence, Chocs; Revue Scientiﬁque et Technique de la Direction des Applications Militaires, CEA, Villeneuve-Saint-Georges, France, 1991.  R.A. Gore, F.H. Harlow, C. Zemach, Computation of mixing layer by a spectral transport model, in: Proceedings of the 8th Symposium on Turbulent Shear Flows, Munich, Germany, 1991, pp. 29-2-1–29-2-6.  D. Besnard, F. Harlow, R. Rauenzahn, C. Zemach, Spectral transport model for variable density turbulent ﬂow, in: Proceedings of the Third International Conference on Compressible Turbulent Mixing, Royaumont, France, 1991.  R.A. Gore, F.H. Harlow, Chemically reacting turbulent mixing layers without heat release, in: Proceedings of the National Fluid Dynamics Conference, Los Angeles, 1992.  D. Besnard, F. Harlow, R. Rauenzahn, C. Zemach, Turbulence transport equations for variable-density turbulence and their relationship to two-ﬁeld models, Los Alamos National Laboratory report LA-12303-MS, 1992.  R. Gore, J. Greene, F. Harlow, Explosive collapse of metal cylinders (U), Los Alamos National Laboratory report LA-12719MS, 1994.  R. Gore, F. Harlow, R. Lohsen, C. Necker, A. Zurek, Explosive collapse of metal spheres (U), Defense Res. Rev. 7 (3) (1997).  M. Steinkamp, T.T. Clark, F.H. Harlow, Stochastic interpenetration of materials, Los Alamos National Laboratory report LA13016-MS, 1995. 432 F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433  E. Scannapieco, F. Harlow, Introduction to ﬁnite diﬀerence methods for numerical ﬂuid dynamics, Los Alamos National Laboratory report LA-12984, 1995.  M. Steinkamp, T. Clark, F. Harlow, Stochastic interpenetration of ﬂuids, I. Model development, Int. J. Multiphase Flow 25 (1999) 599–637.  M. Steinkamp, T. Clark, F. Harlow, Stochastic interpenetration of ﬂuids, II. Numerical solutions and comparisons to experiments, Int. J. Multiphase Flow 25 (1999) 639–682.  R.R. Linn, F.H. Harlow, Wildﬁre Propagation, T-3 Fact Sheet #8, Los Alamos National Laboratory, 1996.  R.N. Silver, R. Gore, J. Greene, F. Harlow, R. Whitman, Radiographic evidence for k-5/3 scaling of density power spectra, Phys. Rev. Lett. 77 (1996) 2471.  R. Linn, J. Reisner, C. Edminster, F. Harlow, J. Winterkamp, W.Smith, J. Bossert, J. Colman, FIRETEC, A Physics-Based Wildﬁre Model, R&D-100 Award, 2003.  T.T. Clark, F. Harlow, R. Moses, Comparison of a spectral turbulence model with experimental data of Rayleigh–Taylor mixing, in: 6th International Workshop on the Physics of Compressible Mixing, in Marseilles, France, 18–21 June 1997.  R.R. Linn, T.T. Clark, F.H. Harlow, L. Turner, Turbulence transport with nonlocal interactions, Los Alamos National Laboratory Report LA-13401-MS, 1998.  R.R. Linn, F.H. Harlow, Transport approach for describing wildﬁres, in: Proceedings of the International Symposium on Computational Technologies for Fluid/Thermal/Chemical Systems with Industrial Applications, ASME/JSME, San Diego, July 1998.  P.N. Wilson, M. Andrews, F.H. Harlow, Spectral measurements of density ﬂuctuations in developing Rayleigh–Taylor mixing, in: 51st Annual Meeting of the American Physical Society Division of Fluid Dynamics (DFD98), Philadelphia, PA, November 22–24, 1998.  P.N. Wilson, M. Andrews, F.H. Harlow, Spectral nonequilibrium in a turbulent mixing layer, Phys. Fluids 11 (1999) 2425.  C.A. Romero, F.H. Harlow, R.M. Rauenzahn, Crenulative turbulence in a converging nonhomogeneous material, Phys. Fluids 11 (1999) 2411.  M.L. Steinzig, F.H. Harlow, Probability distribution function evolution for binary alloy solidiﬁcation, in Proceedings of the Minerals, Metals, Materials Society Annual Meeting, San Diego, California, February 28, 1999.  M.L. Steinzig, F.H. Harlow, Characterization of cast metals with probability distribution functions, in; Proceedings of the Materials Research Society Fall Meeting, Boston, MA, November 30–December 4, 1998.  Michael Steinzig, Francis Harlow, Determination of ﬁnal average grain size in alloy solidiﬁcation, Die Casting Engr. 44 (4) (2000) 48–50.  Duan Z. Zhang, Cheng Liu, Francis H. Harlow, Eﬀects of non-uniform segment deformation on the constitutive relation of polymeric solids, Phys. Rev. E 66 (2002) 051806.  E.N. Harstad, F.H. Harlow, H.L. Schreyer, Single chain stochastic polymer modeling at high strain rates, in: Proceedings of the 2002 ASME Pressure Vessels and Piping Conference, Vancouver, British Columbia, Canada, August 4–9, 2002.  F.H. Harlow, H.J. Shepard, Compilers, theory in action, highlights in the theoretical division at Los Alamos, 1943–2003, Volume I and Volume II, Los Alamos National Laboratory History Report LA-14000-H, April 2003; contains contribution to 14 articles by F.H. Harlow.  M.W. Schraad, F.H. Harlow, A stochastic constitutive law for disordered cellular materials, Int. J. Mech. Sci., to be published. II. Anthropology and Paleontology  J.V. Young, F.H. Harlow, Contemporary Pueblo Indian Pottery Museum of New Mexico Press, Popular Pamphlet Series, 1965 (issued in September 1965).  F.H. Harlow, Tewa Indian Ceremonial Pottery, El Palacio 72 (4) (1965) 13–23.  F.H. Harlow, Historic Pueblo Indian Pottery, The Monitor Press, Los Alamos, New Mexico, 10 February 1967. Reprinted 1968 by Monitor Press. Reprinted 1970 by Museum of New Mexico Press.  P.K. Sutherland, F.H. Harlow, Late Pennsylvanian brachiopods from North-central New Mexico, J. Paleontol. 41 (September) (1967) 1065.  P.K. Sutherland, F.H. Harlow, Pennsylvanian Brachiopods and Biostratigraphy in Southern Sangre de Cristo Mountains, New Mexico, Memoir 27, New Mexico Bureau of Mines and Mineral Resources, Socorro, 1973.  F.H. Harlow, An Enormous Painted Jar from Modern Santo Domingo Pueblo, El Palacio 74 (3) (1967) 44.  F.H. Harlow, The Pottery of San Ildefonso, in: K.M. Chapman (Ed.), The School of American Research, Santa Fe, New Mexico, 1969.  F.H. Harlow, Matte Paint Pottery of the Tewa, Keres and Zuni Pueblos, Museum of New Mexico, Papers in Anthropology, Museum of New Mexico Press, Santa Fe, New Mexico, 1973.  L. Frank, F.H. Harlow, Historic Pottery of The Pueblo Indians 1600–1880, New York Graphic Society, Boston, 1974.  F.H. Harlow, Modern Pueblo Pottery, 1880–1960, Northland Press, Flagstaﬀ, Arizona, 1977.  F.H. Harlow, Glazed Pottery of the Southwest Indians, Am. Ind. Art Mag. 2 (1) (1976) 64–70. F.H. Harlow / Journal of Computational Physics 195 (2004) 414–433 433  F.H. Harlow, Pueblo Indian Pottery traditions, Viltis 36 (3) (1978) 6–7.  F.H. Harlow, An Introduction to Hopi Pottery, half of a booklet published by Museum of Northern Arizona Press, Flagstaﬀ, 1978.  F.H. Harlow, Pueblo Pottery on United States Postage Stamps, The New Mexico Philatelist 31 (2, 162) (1978) 5.  F.H. Harlow, Transitional Pueblo Pottery, 1500–1700, in: Invited Lecture and Distributed Lecture Notes, concourse on Southwest Ceramics, Recursos de Santa Fe, Santa Fe, New Mexico, September 20, 1988.  F.H. Harlow, Rio Grande Glaze Ware, in: Invited Lecture and Published Proceedings, Seminar on Southwestern Native American Ceramics, Millicent Rogers Museum, Santa Fe, New Mexico, 24–26 June 1981.  F.H. Harlow, Pueblo Art: Southwestern Indian Pottery. in: Introduction by Ford Ruthling, The Somesuch Press, Dallas, Texas, 1983.  F.H. Harlow, Two Hundred Years of Pueblo Pottery: The Gallegos Collection, published by Morning Star Gallery, Santa Fe, NM, 1990.  F.H. Harlow, D. Anderson, D. Lammon, The Pottery of Santa Ana Pueblo, School of American Research, Santa Fe, expected publication in 2004.  F.H. Harlow, D. Lammon, The Pottery of Zia Pubelo, School of American Research, Santa Fe, expected publication in 2003.  F.H. Harlow, J. Silverman, Pueblo Indian Pottery, Published by the Silverman Museum, Santa Fe, NM, 2001. (Note: Format is in two forms: books and large folios.).  F.H. Harlow, D. Lammon, The Pottery of Cochiti and Santo Domingo Pueblos, School of American Research, Santa Fe, expected publication in 2005.  Dwight P. Lammon, Francis H. Harlow, Reyes Galvan: A Master Zia Potter, Am. Ind. Art Mag. 28 (2003) 48–59.