Halfway down the Danube

Klytus, I'm bored

God, I hate most current videogames. None of them satisfy my primal need to destroy very large objects in a realistic manner. I spent most of the past week sick, so I decided to use my enforced free time to do something about it.

Turns out that there is enough information available on this Interweb thingy to reverse-engineer the hydrodynamic modeling programs used to simulate Big Things Go Boom physics. You know, shaped charges, nuclear weapons, meteor strikes, the origin of the Moon, et cetera.

It's kind of funny that books on game physics ignore this. There is a niche here to be filled, people!

So let's start large. Let's smash planets.

The simulation of celestial mechanics has been covered extensively in the popular programming literature already. Most often, masses are treated as point sources that "collide" when a minimum distance is reached, in order to prevent crazy numerical divergences. There are sophisticated techniques to get around this problem, for close manuevers and so on. But each planet is tracked individually.

This approach is sometimes called Lagrangian simulation. Each individual mass has its coordinates tracked through space, like a fly on a screen door. But there's another possibility, where each spot on the screen's grid has its mass tracked. This is called the Eulerian method.

Already one can sense a problem with the Eulerian approach: the number of spots on the grid to be analyzed goes up as the size of the volume one is interested in. If you want to see jets of debris extending millions of miles -- and who doesn't? -- you're a little out of luck.

State of the art planet-smashing programs use a modified Lagrangian approach. Planets are composed of bundles of 'kernels', where the mass of each kernel is already spread out over a volume by means of a clever hack. Originally, programs fuzzed the mass using a Gaussian bell curve, the same sort sometimes used to grade exams around an average. But the area under a bell curve is a little time consuming to calculate cleanly, which is why statistics textbooks still have enormous tables in back.

The first simulation of the single impact lunar formation theory -- roughly, another planet strikes the early Earth, and some of the debris consolidates to form the Moon -- published by Benz, Slattery and Cameron in 1986 in the journal Icarus, used a slightly different method. There, the amount of mass in each kernel decayed exponentially from the kernel's center. This was much more mathematically tractable, though the math-savvy will see that a significant fraction of the kernel's mass must extend to infinity. More recent simulations make the math even simpler, using a smooth, quick-to-calculate curve called a spline, that has a sharp cut-off at the kernel's edge.

Another simplification made in the original Icarus paper was to use only the hydrostatic (read: easy) part of the stress tensors in the calculations for the proto-Moon and the proto-Earth. Basically, this means the planets act as giant droplets of liquid. But, as the authors dryly note, "As long as the impact occurs at hypervelocity, shear strength is not important, at least during the impact and for some time afterwards".

The treatment of the shock of impact is interesting. Turns out von Neumann and Richtmyer were looking into this in the 1940s, for what reasons I couldn't begin to imagine. "In the investigation of phenomena arising in the flow of a compressible liquid, it is frequently desirable to solve the equations of fluid motion by stepwise numerical procedures," by which Johnny here means computer programs, "but the work is usually severely complicated by the presence of shocks."

To get around this problem, von Neumann and Richtmyer added a new, non-physical variable to spread the shock out over the sharp edge of impact. "Our idea is to introduce (artificial) dissipative terms into the equations so as to give the shocks a thickness comparable to (but preferably larger than) the spacing of the points of the network." These terms had the same mathematical form as the equations for viscosity, the thickness of a liquid, and so this is sometimes called the artificial viscosity.

Last but not least, we need an equation to describe the behavior of materials under extreme conditions. For some things, the equation is easy. Even in the US, high school kids learn the ideal gas equation:

P * V = n * R * T

where P is the pressure of a gas, V the volume, n the number of moles (don't ask), R a constant, and T the absolute temperature. This can be tweaked to:

P = (gamma - 1) * rho * E

where rho is the density of the gas, E the internal energy, and gamma for an ideal gas exactly 5/3. It turns out you can use this for the Sun and only be a percent or two off.

Of course, planets, even when they are simulated as giant droplets of liquid, don't behave like ideal gases. At least, of course, until they are vaporized by impact. Turns out a guy named Tillotson in 1962 hacked together a reasonably fast formula to deal with regular Joe solids that might suddenly vaporize. It was first published in a General Atomics report, I couldn't say why.

Tillotson's formula comes in two parts, the first for when the internal energy of the material is less than the energy needed to begin vaporization (the cold part), and the second after the internal energy has passed the point of complete vaporization (the hot part).

The cold part is:

P = a * rho * E + (b * rho * E / (E / (E_0 * eta ^ 2 + 1))) + A * mu + B * mu ^ 2

Apologies for the alphabet soup. P, rho, and E are as before; eta is the current density divided by the standard density; and mu is that number minus 1. E_0, a, b, A, and B are all material dependent constants.

The hot part is somewhat worse:

P = a * rho * E + ((b * rho * E / (E / (E_0 * eta ^ 2 + 1))) + A * mu * exp(-beta * ((1 / eta) - 1)) * exp(-alpha * ((1 / eta) - 1) ^ 2)

where alpha and beta are also material dependent constants that control how quickly the material begins to act like an ideal gas.

Pressures in between are extrapolated linearly, easy enough.

Note the parts I have put in bold. They're similar to the above equation for an ideal gas, no? Theoretically, a should be exactly 2/3. Observationally, for most materials, a = 0.5 works better. Also, some tweaking shows that A behaves very similar to a bulk modulus, so you can pull out a speed of sound = sqrt(A / the standard density).

Still, most of these coefficients are a little difficult to calculate (or frankly, even to find). So I've compiled a list, taken mainly from Melosh's Impact Cratering. The current state of the art goes straight to the thermodynamic state variables, but perhaps some game designer or recreational programmer might find these of use.

For aluminum, the standard density is 2700 kg per cubic meter, A 75.2 gigapascals, B 65 GPa, E_0 5 megajoules per kg, the energy of incipient vaporization 3.0 MJ/kg, the energy of complete vaporization 13.9 MJ/kg; a 0.5, b 1.63, alpha 5, and beta 5.

For basalt, the standard density is 2700 kg/m^3, A 26.7 GPa, B 26.7 GPa, E_0 487 MJ/kg, the energy of incipient vaporization 4.72 MJ/kg, the energy of complete vaporization 18.2 MJ/kg; a 0.5, b 1.5, alpha 5, and beta 5.

For granite, the standard density is 2680 kg/m^3, A 18 GPa, B 18 GPa, E_0 16 MJ/kg, the energy of incipient vaporization 3.5 MJ/kg, the energy of complete vaporization 18 MJ/kg; a 0.5, b 1.3, alpha 5, and beta 5.

For iron, the standard density is 7800 kg/m^3, A 128 GPa, B 105 GPa, E_0 9.5 MJ/kg, the energy of incipient vaporization 2.4 MJ/kg, the energy of complete vaporization 8.67 MJ/kg; a 0.5, b 1.5, alpha 5, and beta 5.

For limestone, the standard density is 2700 kg/m^3, A 40 GPa, B 67 GPa, E_0 10 MJ/kg, the energy of incipient vaporization 2.5 MJ/kg, the energy of complete vaporization 14 MJ/kg; a 0.5, b 0.6, alpha 5, and beta 5.

For water ice, the standard density is 917 kg/m^3, A 9.47 GPa, B 9.47 GPa, E_0 10 MJ/kg, the energy of incipient vaporization 0.773 MJ/kg, the energy of complete vaporization 3.04 MJ/kg; a 0.3, b 0.1, alpha 10, and beta 5.

(The B value for water ice is extrapolated; it possibly should be significantly higher. Also, note that a is low. This is due to the effects of molecular water. You'd expect this effect in ammonia or methane ice as well.)

Alas, I don't have any values on hand for exotic materials like diamond or metallic hydrogen (the latter, I am not sure if anyone does), so smashing open Jupiter is a little difficult.

But smashing terrestrial planets open, or even colliding stars, should be well within the reach of current videogame platforms. L.B. Lucy in 1977 simulated a protostar spinning apart with 300 particles in four minutes on an IBM 360. Benz, Slattery and Cameron used 2048 particles in their lunar impact simulation less than a decade later. Surely Moore's Law has caught up to the point where this could be done on an Xbox or a Play Station?

Posted by coyu at December 26, 2004 03:06 AM