SPAM-BASED RECIPES FOR CONTINUUM SIMULATIONS

Smooth Particle Applied Mechanics (SPAM) is a technique which provides a versatile approach to simulating many difficult problems in continuum mechanics, such as the breakup of a cavitating fluid and the penetration of one solid by another. SPAM also provides a simple evaluation method for all the continuum variables, as well as the spatial gradients required by the evolution equations. This global evaluation simplifies interpolation, reasoning, and Fourier transformation. Because SPAM is so flexible and easy to program, it should be included in the toolkit of anyone doing simulations.

assumed to vary continuously through-uids, and solids, both near to and far out the space-time region of interest. from equilibrium. Flows in systems as Undergraduate texts describe the solu-diverse as phase-separating fluid mixtion of "easy" continuum problems, in-tures, deforming metals, breaking cluding the flow of heat in response to rocks, and colliding astrophysical bodimposed thermal boundary conditions ies can all be treated with SPAM. and the linear-elastic response of solids Although the dependent variables p , to imposed loads. A variety of grid- 3 , and e are defined throughout space based approaches work very well for and for a continuously varying time, such problems. For more difficult prob-the solutions of the partial differential lems, in which the structure under study equations (Equations 1 through 3) can undergoes extreme shape changes and be made tractable by converting them forms new surfaces while resisting in-p dt p (') to ordinary differential equations, terpenetration at interfaces, more so-spanning space with a grid of points, and computing the time histories for a phisticated and flexible techniques are p-p.' = -v. P required. dt closely spaced set of discrete times. A is such a technique.'-6 SPAM provides p-= dt @=-p:V3-V'Q. (3) best choice for the time integration. a versatile approach to simulating Until Leon Lucy' and Joseph Monmany difficult problems in continuum Vectors are denoted by an arrow and aghan' conceived of SPAM in 1977, mechanics, such as the breakup of a tensors are in boldface. If we denote spatial gradients of the field variables cavitating fluid and the penetration of the components of a vector as A>,then p, 3 , and e were typically evaluated on one solid by another. SPAM also pro-the components of the tensor AB are an underlying spatial grid of interpolavides a simple evaluation method for tion points. (Joseph Monaghan of all the continuum variables, as well as Monash University, Adelaide, Austhe spatial gradients required by the tralia, is a prolific and creative expoevolution equations. This global eval-nent of smooth-particle methods. The uation simplifies interpolation, rezon-reader is encouraged to seek out Moning, and Fourier transformation. Be-aghan's recent work on the Internet.) cause SPAM is so flexible and easy to Both comoving Lagrangian grids, program, it should be included in the moving with the material, and stationtoolkit of anyone doing simulations.
ary Eulerian grids, fixed in space, were

SPAM
The primary application of SPAM is to the solution of problems in continuum mechanics, where the governing partial differential equations describe the evolution of the mas density p, velocity v' , and energy per unit mass e in terms of gradients of the velocity, pzessure tensor P, and heat-fluxvector Q : Smooth Particle Applied Mechanics de used in this way. In both cases, the spatial derivatives appearing in the continuum equations were approximated as simple finite differences of field-point values at neighboring points. Although this interpolation approach works well for simple problems, it has difficulties with the chaotic irregular flows typified by mixing, penetration, fracture, and turbulence. The distortion of a comoving mesh, the difficulty of following material interfaces, the prevention of overlaps boundary conditions all can cause major headaches. SPAM solves the problem of choosing a spatial grid by introducing smooth particles, with a finite spatial extent or the set of points used to interpolate the field variables' values. The interpolation grid varies continuously as the particles move. This spatial interpolation method is in itself a useful algorithm for smoothing data known on an irregular set of points. To illustrate the interpolation method, consider a 2D space with a set of points (the particle coordinates), all of which have the same mass m and the same weight function w(r). The particle mass density (or probability distribution) can be thought of as spread out in space according to w(r) with a range h characterizing the size of es. The sum of all particle contributions to the p at a point 7 gives the total density (or prob-which are defined everywhere, not just at the particle positions 5. In particular,& differs from f( 7 = <) if, as is usual, there are local fluctuations in$ For reasonable accuracy at manageable cost, the range of the weight function w should be such that the weighted averages include contributions from about 20 particles. Lucy introduced a useful form for ~( r ) : ' within a fixed mesh, and the smooth implementation of (6) width. The instantaneous coordinates of the particles define is chosen to make the representscase of an infinite number of interpolating particles.
For r just inside the range h, the weight function in Equation 6 vanishes as the cube of the distance to the cutoff, m(r) = (hr)j. Thus w has two continuous derivatives everywhere. Likewise, the first and second spatial derivatives of field-variable sums such as ability). We define the mass density as have no discontinuities. Such a smooth-particle i tion method can be applied to molecular dynamics, where the individual particles have kinetic and potential energies. In this case, smooth-particle interpolation provides a twice-differentiable, continuous representation of the two-particle energy functions based on the summed contributions from discrete points. This representation is much more useful than the usual delta-

Continuum equations
The SPAM representation, with two continuous derivatives for the field variables, is the key to solving the continuum equations. Consider P,, the total mass density at the location of particle i. It is given by a special case of the general rulefor p(?),with ?+<: where the sum is over all particles within the range h of 7 and m is the particle mass. Note that the form of Equation 4 is similar to the usual expression for the density, but with the usual delta functions replaced by a "smooth" (at least twice differentiable) weight function w. This is the reason the particles are referred to as smooth particles. For a gen-era1 field variable5 the same procedure of adding all the tributions motivates the fundamental f(r')P(?) = m.C&-r,i). At any position r' , the interpolated value of the field variable f(7) is given by a weighted average of contributions& from nearby particlesj within the range h. We emphasize that the discrete particle properties$, one for each particle, differ from the continuum of interpolated averages f(7), Note that the choice of the normalization of w guarantees that the integral over all space, JJp(x,y) dx dy, equals the total mass mN, where Nis the number of smooth particles.
Not only mass but also momentum and energy are conserved by SPAM. The latter two conditions follow directly from the smooth-particle version of the continuum evolution equations for velocity and energy. To see this, we use the definition (Equation 5 ) to find the smooth-particle representation of spatial gradients: 3 Because the4 are particle properties rather than spatial averages, the gradient operator affects only the weight function w through its explicit dependence on r' . f = (P/p2), we can show that the smoothparticle equations of motion conserve momentum exactly. We first evaluate the divergence of (P/p) in the usual con-By setting f ng the rules of ordinary calculus:

V . ( P / p ) = -( P / p ' ) . V p + ( l / p ) V . P .
(10) It is important to note-as will be emphasized in the following example-that the special case of a scalar pressure P proportional to p2 gives equations of motion isomorphic to those of molecular dynamics, with the weight function w playing the role of a pair potential. Whenever the quotient (P/$) varies slowly in space, the smooth-particle continuum dynamics resembles ordinary molecular dynamics. Even in the most general case of a tensor pressure, the smooth-particle equations of motion conserve momentum exactly. To see this, consider the time variation of the sum of the particle momenta, Because the contributions of each interacting { i ,~] pair exactly cancel in this sum, the sum must vanish. Thus the smooth-particle equation of motion (Equation 13) conserves momentum. The reader should be able to verify, using a similar argument, that the continuum energy equations can be written in a smooth-particle form that conserves energy exactly: If we combine the two gradients and replace F by < , we find the smooth-particle form for the acceleration at the POsition of particle i. This form is the equation of motion for particle i: 4. Compute P, and (z, at each particle, using the constitutive relations, which include the equation of state P@, e) or P@) as well as the phenomenological relations governing diffusive, viscous, and conductive flows.

5.
Compute G, and e, using Equations 13 and 14, including the pressure tensors and heat flux vectors computed in step 4.
integrating the differential equations for g / d t =   The corresponding expressions for the velocity and temperature gradients at the position of particle i guarantee that the gradient contributions of each pair of particles are proportional to the corresponding velocity and temperature differences: To illustrate the smooth-particle method we apply it to a relatively difficult problem, the free expansion of a 2D ideal gas, an Euler fluid with vanishing transport coefficients, into a container four times its original s i~e .~,~ The equilibrium proportionality constant a function of energy only. A conventional fixed-mesh Eulerian grid approach to this problem exhibits catastrophic instabilities unless artificial damping is added. A conventional Lagrangian approach, with its mesh elements moving with the fluid, exhibits two kinds of numerical shear instabilities, the butte$y and the hourglass, as shown in Figure 1. Both instabilities can be troducing numerical viscous damping terms to but the programming effort is severe. The sm method requires no unusual precautions. No catastrophic instabilities can occur, because the particle grid is continually evolving, with the field variables everywhere smooth. no transport coefficients to increase the entropy above its initial value, it is a challenge to understand the entropy increase that must result, AS = Nk In 4. The resolution of this problem arises naturally from the SPAM simulations?16 For the free-expansion problem, it is straightforward to choose the initial particle coordinates (a square grid is fine). By introducing small random offsets-equivalent to small initial random particle velocities-the square symmetry of the grid can be broken. In the simulations, we choose the smooth-particle mass and initial density both equal to unity and use the isentropic equation of state heat flux vector 6 vanishes. Thus the smooth-particle representation of a special ideal gas isentrope gives particle trajectories identical to those of molecular dynamics. The reader should be able to show that the continuity equation (Equation 1) and the energy equation (Equation 14) both lead to the same result: Because these two tinuum equations are simultaneously satisfied by the density definition (Equation 4), we do not need to integrate the energy equation. Figure 2 shows snapshots of the particle motion during the free-expansion simulation. The underlying equations of motion have simple cubic forces derived from the weight function given in JQpt i on 6. The range h of the weight function w was chosen to be equal to six times the initial square-lattice nearest-neighbor spacing of unity. The density contours-equivalent to internal energy contours-and the kinetic energy contours shown in Figure 2 were computed with the smooth-particle definitions. The smooth-particle approach captures velocity fluc-4 = 4 = x (~~ -~~).v,~(17; -50 . io The velocity fluctuations are crucial to an understanding of the entropy increase in &IS problem. In equilibrium s, the entropy density (per unit volume) for a 2D ideal gas is equal to -k(p/m) In (dp), where k is Boltzmann's constant. In the nonequlibrium case considered here, the local velocity fluctuations, relative to the mean flow, must be included in the local internal energy. The smooth-particle approach gives these fluctuations-absent in the usual grid-based approaches-at every point. The nonequilibrium entropy density, including this local kinetic energy density due to fluctuations, isSs6 expected increase of the entropy as predicted from statistical mechanics, AS = Nk In 4.
The results shown in Figures 2 and 3 were computed with periodic boundary conditions. Mirror boundary conditions, applied as illustrated in Figure 4, are another possibility and lead to very similar results. With mirror boundaries, each particle near, but inside, the boundary (within the range h) produces an image particle outside the boundary, to which boundary values of velocity and energy per unit mass can be attributed. The treatment of problems in which the location of the boundary is unhown (as in the surf near a beach or the formation of a crack or tear) is a challenging research area for smooth particles. Now that parallel processing with thousands of processors is a reality, some of the difficulties in treating boundaries are being alleviated by more computer power.
Although the example discussed here involves identical particles with equal masses and the same weight function, there are circumstances that justify the extra work of a more general treatment.* It is relatively easy to use a positiondependent particle size, with h depending upon the local density. A more sophisticated treatment introduces ellip- -k(p/m)ln[(e+3((V-(v))2))/p].

(22)
The total entropy for the smooth particles can be evaluated in either of two ways: by integrating this expression over a regular grid spanning the system, or by summing up the individual particle contributions. The two methods agree quite well (see Figure 3). In only a little more than the time required for sound to travel across the system, both the resulting integral and the corresponding sum reproduce the ---soidal particles, with the principal-axis lengths depending on the density gradient.
An underlying set of smooth particles makes it possible to evaluate field variables on a regular grid, not just at the particle locations. This flexibility is particularly advantageous if a rectangular mesh is required, as is the case for some graphics programs and for computing fast Fourier transforms of the field variables.
Rezoning can be accomplished easily with smooth particles. In the event that more detail is required in a particular region, a particle can be replaced by two or more smaller ones, choosing the new masses, velocities, and energies so as to conserve mass, momentum, and energy. Similarly, if a region has too many particles, two can be combined into a more massive and energetic single particle. The method can be extended in many ways to treat special situations. Chemical reactions can be introduced. Electromagnetic fields can be included by using tree techniques to evaluate the effect of long-range forces. As is usually the case, the major advances in numerical simulation methods have resulted from the desire t o simulate challenging problems. : :