CS205b/CME306Lecture 31 Smoothed Particle Hydrodynamics (SPH)1.1 Representation and SimulationIn the previous lecture, we showed how to define density throughout space by smearing out massesat discrete points in space.ρ(x) =XimiW (x − xi)We did not describe how one might choose those m asses mior the locations xiwhere they live.If we instead have some initial conditions, we would like to choose a suitable set of miand xisothat the dens ity profile approximates the initial conditions. Given a certain number k of particles,we can choose k values miand k values (in 1D) xigiving us 2k variables. We can write down2k equations by m aking 2k measurements of ρ, and solve for a suitable set of initial masses andlocations, though run into issues with overdetermined/underdetermined systems.This gives us some connection between the continuous representation and a discrete r ep resen-tation. The idea of numerical simulation is to evolve the discrete version forward in time. Thecontinu ou s one is what we see in the real wor ld, and it evolves f orward in time based on somephysical laws. If at some later time we could compare the the continuous one to the discrete one,we would like to have the two match in some way. There will be errors between the two, whicharise from two sources: errors in the initial discretization and errors that accumulate during thecourse of the simulation. In this way, we may use numerical simulation to pr ed ict the outcome ofa continuous system.1.2 Other AttributesWe are not, however, limited to attaching density to our chunks. We can attach any attribute Aito the chunks 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 volume weighting that cancels out the weighting in W (x). If we integrate thecontribution of one chunk throughout space,Z∞−∞A(x) dx =Z∞−∞miρiAiW (x − xi) dx =miρiAiwhich is just the volume weighted attribute as one would expect.1Because 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 attributes around. We automatically conserve mass with this method, so the nextstep would is to consider momentum.2 Conservation of Momentum (Forces)At this point we would like to actually simulate something. We have not yet introduced any meansby which this motion can be computed, and this is where forces are introduced.Newton’s second law p rovides the neces s ary relationship between forces and motion and may bewritten as F = ma = p′, where p is the momentum of a particle. One may also view this relationshipas an extension to conservation of momentum. Note that Newton’s first law is conservation ofmomentumPipi= constant, which a consequence ofPiFi= 0, a system experiencing no externalforces. Newton’s third law (equal an d opposite reactions) requ ir es that f orces between particlesoccur in equal and opposite p airs, so in the absence of external forces, net force is still zero and themomentum of the s ys tem is conserved. Newton’s third law provides a convenient means for enforcingconservation of momentum in the Lagrangian framework (we will use this when formulating springslater).Particles may now be evolved through space as long as the forces acting on them can becomputed. One of the simp lest and most important forces is gravity. For our purposes, gravitymay assumed to be constant throughout space. Then, we compute the force on a particle dueto gravity as F = mg, so that a particle experiencing no other forces simply falls with constantacceler ation a = g. A somewhat more interesting force is a simple drag force F = −kv, where k isconstant and v is the particle’s velocity. A particle with higher velocity feels more drag, and theresulting force opposes its motion. A particle experiencing only this force slows down but neverreaches rest or changes its direction.3 Linearized SystemForce is in general a function of both the positions and velocities of the particles in a system.That is, F = F (x, v). In the interests of writing down a linear system to analyse stability, it isconvenient to approximate this force as F (x, v) ≈ F (x0, v0) + Fx(x, v)x + Fv(x, v)v and also ignorethe inhomogeneous term F (x0, v0), since Duh amel’s pr inciple states that the inhomogenous termdoes not affect stability. Note that this approximation omits gravity, since it is inhomogeneous.We can typically make these simplifications when looking at stability. Forces of the form Fx(x, v)xare r ather like spring forces, and forces of the form Fv(x, v)v behave like damping forces. UsingFx= maxand Fv= mav, the motion of the particle is described by the fi rst order linear system xv!′= 0 1axav! xv!.2The eigenvalues of this system areλ =av±pa2v+ 4ax2and have units of Hertz (s−1). Solutions look eλt, so well posedness requires the real part of theeigenvalues be nonpositive, Re(λ) ≤ 0. This places some restrictions on the way we can modelforces of nature to prevent the system from blowing up.It is necessary for ax≤ 0. If a particle experiences n o force at the origin but experiences astronger force as it moves away, that f orce 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 period behavior. If av< 0, the system has exponentialdamping. If av= 0, the system is undamped. Note that in the undamped case, the eigenvalues arepure imaginary. In the underdamped case |λ| =√−axdoes not depend on the amount of dampingapplied.When −av> 2√−ax, we call the system overdamped. The eigenvalues are real and distinct,and the solution exhibits exponential decay only.When −av= 2√−ax, we call the system critically damped. The eigenvalues are equal, and thesolution exhibits exponential decay with at m ost one