CS205b/CME306Lecture 21 Smoothed Particle HydrodynamicsIf we consider mass as attached to chunks, we automatically conserve mass as we move th e chunksaround. However, because the laws governing motion often involve derivatives of quantities storedon the chunks, it is us eful to have a definition of these properties everywhere in space that issufficiently differentiable. This can be achieved by spreading attributes associated with chunks,such as their mass, over some local region of space. This is the basic idea behind the SPH method.Let W (x) be some sufficiently differentiable function that can be used to spread attributesover space. In p articular, consider a chunk centered at xiwith mass mi. Then define the densitycontribution from this chunk at any point x in space to be ρ(x) = miW (x −xi). Because we wouldlike to conserve mass, we insist that the total mass spread throughout space is the m ass assignedto that chunkmi=Z∞−∞ρ(x) dx =Z∞−∞miW (x − xi) dxwhich leads to the requirementZ∞−∞W (x) dx = 1.It is also desirable for W (x) to be symmetric about the origin (so that dens ity is spread evenly inall directions about its center of mass). It is also more efficient if W (x) has local influence in thatW (x) = 0 everywhere outsid e some region around the origin. Note that W (x) has units of one overvolume, since it yields unity when integrated over the volume of space.To be very useful, our model must contain many chunks. In this case, we d efi ne the den s ity atany point as the sum of the density contributions of every chunk as illustrated in Figure 1.ρ(x) =XimiW (x − xi)bx1bx2bx3Figure 1: Three chunks located at x1, x2, and x3. The black curves show the density contributionof a single chunk at each point in space. The dotted line shows the density profile of space.1We are not limited to attaching density to our chunks. We can attach any attribute Aito thechunks and use W (x) to distribute it. This is typically done usingA(x) =XimiρiAiW (x −xi)where ρi= ρ(xi). Note that we get back our definition of ρ(x) if we let Ai= ρi. The extrascaling factor is a volu me weightin g that cancels out the weighting in W (x). If we integrate thecontribution of one chunk throughout s pace,Z∞−∞A(x) dx =Z∞−∞miρiAiW (x −xi) dx =miρiAiwhich is just the volume weighted attribute as one would expect.Because we have a smooth definition of A(x) everywhere, we may compute its derivatives. Forexample, ∇A can be computed as∇A =XimiρiAi∇W (x − xi).We have now described how to represent mass, other scalar quantities, and derivatives of thesescalar quantities in space just based on the idea of having quantities attached to attributes andmoving these attr ibutes aroun d. We automatically conserve mass with this method, so the nextstep would is to consider momentum.2 ForcesUp to this point we described how chunks of material evolve in time and how the density of anarbitrary point in space could be computed u nder a scheme such as Smoothed Particle Hydrody-namics model. In this way, any property could be defined and differentiated over space. However,we have assum ed u p to this point that the m otion of p articles thr ough sp ace is known. We havenot yet introduced any means by which this motion can be computed, and this is where forces areneeded.Newton’s second law p rovides the n ecessary relationship between forces and motion and maybe written as F = ma = p′, where p is the momentum of a particle. One may also view thisrelationship as an extension to conservation of momentum. Note that conservation of momentumPipi= constant is a consequence ofPiFi= 0, a system experiencing no external forces. Notethat Newton’s third law (equal and opposite reactions) requires that forces between particles occurin equal and opposite pairs, so in the absence of external forces, net force is still zero and themomentum of the system is conserved. Newton’s third law provides a convenient means for enforcingconservation of momentum in the Lagrangian framework.Particles may now be evolved through space as long as the forces acting on them can becomputed. One of the simplest and most important forces is gravity. For our purposes, gravitymay assumed to be constant throu ghou t space. Then, we compute the force on a particle dueto gravity as F = mg, so that a p article experiencing no other forces simply falls with constantacceler ation a = g. A som ewhat more interesting force is a simple drag force F = −kv, where k isconstant and v is the particle’s velocity. A particle with h igher velocity feels more dr ag, and theresulting force opposes its motion. A particle experiencing only th is force slows down but neverreaches rest or changes its direction.23 Linearized SystemForce is in general a function of both the positions and velocities of the particles in a system. Thatis, F = F (x, v). In the interests of writing down a linear system, it is convenient to approximatethis force as F (x, v) ≈ F (x0, v0) + Fx(x, v)x + Fv(x, v)v and also ignore the inhomogeneous termF (x0, v0). Note th at this approximation omits gravity, since it is inhomogeneous. We can typicallymake these simplifications when looking at stability. Forces like Fx(x, v)x are rather like sprin gforces that get s tr onger as the , and forces like Fv(x, v)v behave like damping forces. Using Fx=maxand Fv= mav, the motion of the particle is described by the first order linear system xv!′= 0 1axav! xv!.The eigenvalues of this system areλ =av±pa2v+ 4ax2and have units of Hz = s−1. Solutions look eλt, so well posedness requires ℜλ ≤ 0. This placessome restrictions on the way we can model forces of nature to prevent the system from blowing up.It is necessary for ax≤ 0. If a particle experiences no force at the origin but experiences astronger force as it moves away, that force should be a restoring force rather than one that push itaway harder as it moves farther away and causes exponential growth.Similarly, it is necessary for av≤ 0. A force that did not satisfy this would tend to apply forcesin the direction of motion that get stronger as the particle moves faster and result in exponentialgrowth.When −av< 2√−ax, we call the system underdamped. The eigenvalues contain imaginarycomponents, and the solution exhibits perio d behavior. If av< 0, the sys tem has exponentialdamping. If av= 0, the system is undamped. Note that in the undamped case, the eigenvalues arepure imaginary.When