Can bullets fired into the air kill?
Since joining the military I'd always wondered about this—can a bullet fired into the air, be it intentional or not, drop back down with enough speed to kill? During my training we talked about whether it was possible, and sometimes we would postulate that it was, and sometimes not; and in the first place those were conclusions drawn from the imaginations of bored men, tainted with as much uncertainty as we had at our own future in the army.
This is an attempt to satisfy my curiosity—with some physics and some simulations, and some graphs to look at.
It will be an absurdly long post. I'll provide a rough road-map right here, so that navigating this hunk would be easier for you. We will begin by revisiting the fundamental physics of projectile motion, with and without air resistance. Whilst doing so there'll be several digressions on interesting observations, which I hope would provide a more comprehensive theoretical coverage of projectile motion.
Following that, we'll introduce several drag models, obtained from literature, ballistic calculators, and reference charts. These would all be in the form of velocity-dependent drag coefficients, except for a single constant-coefficient model that we keep for the sake of comparison. We'll also run a finite-element computational fluid dynamics simulation ourselves, and obtain a drag model of our own.
To decide if the bullet would indeed be lethal, we evaluate the velocities of the projectiles under various impact conditions, and also read up on some terminal ballistics literature to gauge their wounding ability. By this point we would have answered our question., but we bring it a step further. We'll proceed to analyze the probability of a misfired bullet lethally hitting someone, which would involve computational geometry and some data mining, in addition to the physics of bullets.
I'll kill the suspense now—the conclusion seems to be that in most situations the bullet would not kill, for it slows down due to air drag much too fast. Under certain conditions there is a chance of lethality, but in an overwhelming majority of cases, the bullet is harmless. At the same time, we also find that the probability of the bullet hitting someone is astronomically low. And so firing a bullet blindly into the sky might not be so catastrophic after all.
Projectile Motion and Interesting Things
Without air resistance, a projectile launched at some angle θ from the horizontal would follow a parabolic trajectory. The equation of motion of the projectile can be written as follows, where we have adopted a Cartesian coordinate system with the projectile launched at the origin. Initial conditions are also written.
This set of equations have an analytical solution, of course.
This is the familiar parabolic projectile motion as taught in classrooms. We know that this is merely an idealization, for air resistance is present in the real world. For bullets travelling at several times the speed of sound, the drag force will be enormous—and so it is paramount that we take it into account. The retarding effect of air resistance is the reason why rifle projectiles have an effective range of at most a few kilometers.
We can consider Stokes' drag, which is linear in velocity. This linear drag can be used to characterize the viscosity of liquids, and is also frequently invoked to analyze the swimming behaviour of small aquatic organisms, or to understand sedimentation under gravity or in a centrifuge. It is indeed quite useful in the real world.
Above b is the drag coefficient, which depends on the geometry of the object and the properties of the fluid. The initial conditions are the same for parabolic motion, with specified launch velocity and launch angle. This set of equations, too, have analytical solutions.
The trajectory is, of course, no longer a parabola; we have additional exponential factors. We will see graphically later that this set of solutions describe a somewhat distorted trajectory, with the downward leg being steeper than the upward leg.
And there's also Newton drag, which is quadratic in velocity.
The homogeneous form of these equations—corresponding to a zero-gravity scenario—does have an analytical solution; but the non-homogeneous form does not. This is not much of an issue, because we can proceed via numerical means. Mathematica, a mathematical symbolic computation program, does allow us to solve these equations using a fourth-order Runge-Kutta numerical scheme pretty easily. And so we can evaluate the trajectories of projectiles and plot them out for visualization.
We notice that quadratic drag has a much larger effect on projectiles of a given weight, as compared to Stokes' drag. And we can also observe the asymmetry in the trajectories, as introduced by air resistance.
Digression #1: Flight Time of Projectiles
One thing which is not immediately obvious in the graph above, is that projectiles under the influence of air drag actually takes less time to reach the ground again, as compared to those without air drag. It is an interesting observation, because we'd have intuitively expected air drag to slow them down—which it indeed does—but the vastly shorter range caused by air resistance far outweighs the retardation. I'd generated an animation, to show this phenomenon. Pay attention to the landing times of the projectiles.
Is this true in general? Would projectiles under the influence of air resistance always land earlier than those without drag, regardless of launch velocity and launch angle, and projectile mass and drag coefficient? We can convince ourselves, somewhat, by generating a graph of flight time—the time elapsed before landing—against launch velocity, and animate it for different launch angles and drag coefficients. These animations are presented below. We can see clearly in these animated graphs that projectiles with air resistance land earlier; and this seems to be true for all cases.
So it seems to be true in general. Now the next logical question to ask—is there a way to prove this, mathematically? And there is, indeed, at least for the case of Stokes' drag.
We will begin by finding the flight time of a drag-less projectile, which is quite trivial. This will be the baseline for comparison.
We can also find the flight time for a projectile under the influence of Stokes' drag, albeit with a little more work.
Above we have introduced the Lambert W function, denoted W. What the Lambert W function is, and its inherent properties, are a little complicated, and we would not touch on these details here. Our aim is to prove that the flight time of the Stokes' drag projectile is less than that of the drag-less projectile, and we can perform some manipulation right away to simplify this postulated inequality.
Now we make the observation that equality is achieved when z = 1.
And therefore, if for all z from z = 1 onwards the W function does not increase more quickly than (z - 2), then we would have proven our statement. In other words, we want the derivative of W to be smaller, or equals to, unity.
This is equivalent to the following condition, which we want to show is true.
And this is necessarily true, because z cannot be smaller than unity by our definition. It is helpful to note that W is a monotonously increasing function, which allows us to do the arithmetic as presented.
And this completes the proof. The mathematics does indeed show that a projectile under the influence of Stokes' drag would land before a projectile without air drag, given that they are launched at the same angle and same initial velocity; and this remains true for all launch angles and velocities.
We would not attempt a proof for the quadratic drag case, because numerical solutions and direct mathematical analyses don't go well together. There is, however, a paper titled "Approximate Analytical Description of the Projectile Motion with a Quadratic Drag Force", by Peter Chudinov; and in it he presented some approximate analytic solutions for the quadratic drag projectile. In particular, we can extract his results for projectile flight time, and compare it to the drag-less case.
This approximate analytic solution by Chudinov is quite accurate, with typical errors of 1-2%. Even so, the fact that it is an approximation means that this is by no means a valid proof. Notably, there are exact solutions for velocity expressed as a ratio of infinite series, as published in "An analytic solution to the equations of the motion of a pointmass with quadratic resistance and generalizations" by Ray and Fröhlich; but I had been facing difficulties in using their solutions to prove the flight time inequality.
I had included this digression because I felt this whole air-resistance-makes-projectiles-land-faster thing goes very much against our intuitions, and it deserves some attention. And, quite interestingly, the use of Lambert W functions in projectile motion seems to be a fairly recent development, first published in 2004 by Packel and Yuen, in their paper "Projectile Motion with Resistance and the Lambert W Function". They said in their paper that it wasn't their idea—it was Mathematica's—which is quite awesome. I'm really excited to see what impact symbolic computation would have on the research scene of the future.
Digression #2: Optimal Launch Angle
Here's another digression. Again, I had chosen to include this because I think it's interesting. Our theoretical discussion on projectile motion just wouldn't seem complete without it. The actual analysis on whether or not bullets shot into the air can kill will resume in the next.
We know that the optimal launch angle for a projectile without air resistance is always 45 degrees to the horizontal. The derivation for this is presented, for a quick recap. Below we use R to represent projectile range.
It is noteworthy that the optimal angle is a constant—that it does not depend on launch velocity, or the mass of the projectile. We expect this by intuition, but this is a special case. With air resistance, the optimal angle ceases to be a constant; it will be dependent on launch velocity, projectile mass, drag coefficient, etcetera—and this dependency will not be a simple one.
We'll take a look at the Stokes' drag case. Writing out the equations in full would certainly be unpleasant, so we had made use of z as introduced earlier.
And of course attempting to solve for optimal launch angle analytically would lead to certain death, so we make use of the good old Newton-Raphson numerical method instead. We plug in projectile mass on the order of a few grams, which is typical for bullets. We can then plot out the optimal launch angle for different drag coefficients, as shown below.
And so we see that with air resistance, the optimal launch angle for maximum range is always below 45 degrees; and this angle decreases as air resistance gets stronger. Intuitively, we can reason that the path length traveled by the projectile gets greater as its launch angle increases, and so the retardation brought about by air resistance will increase in extent, thus resulting in a decrease in range. And so, to have maximal range, we have to bring the launch angle down.
A similar approach can also be taken for the quadratic drag case, but because the equations of motions cannot be solved analytically, we'd have to employ numerical differentiation from the start. The profile of optimal launch angle is noticeably different from that of Stokes' drag.
And again we see that the optimal angle is always below 45 degrees—and the same intuitive reasoning can be applied to make sense of this. This short digression does well to illustrate the complexity and intricacies involved in external ballistics. It is not surprising that scientists and militaries throughout the ages had come up with ingenious solutions to calculate projectile trajectories, first giving life to tediously compiled ballistic charts, and then to mechanical calculators, and now to complicated pieces of software.
That's all for this digression. We'll move on to analyzing the bullets.
Reynolds Number and Quadratic Drag
We mentioned Stokes' drag and quadratic drag, and we did a fair amount of analysis on each; but in reality only quadratic drag is useful for our purposes. This is because Stokes' drag is an appropriate model for low-velocity flows, and clearly, bullets are not low-velocity objects. The way we differentiate between low- and high- velocity flows is through an approximate metric called the Reynolds number, the ratio of viscous to inertial forces in the flow.
Above ρ is the density of the fluid, v is the fluid velocity, L the characteristic length, and μ the dynamic viscosity of the fluid. Stokes' drag applies only when the Reynolds number is significantly less than unity. In our case, an approximate calculation at a bullet velocity of 800 m/s yields a humongous Reynolds number. Do note that for most rifles, the muzzle velocity is in excess of 800 m/s.
Even for a bullet that has slowed considerably due to air drag, the associated Reynolds number would still be much above unity. This dictates that we should be using quadratic drag. Furthermore, it indicates that there is likely to be turbulence in the flow, especially in the trailing regions; we'll talk more about this later.
We should also note that while the drag coefficient is usually understood to be a constant—hence the name coefficient—it really is not. There's a tonne of factors that affects how drag scales with velocity, amongst them the points of flow separation, turbulent effects, transonic transitions, and inherent non-linearity in the fluid. A constant coefficient might be a reasonable simplification for a narrow range of subsonic velocities; but because we're dealing with a wide range of velocities across both subsonic and supersonic regimes, we have to consider the actual non-constant coefficient functions. This shall be the focus of the section following next, "Bullet and Drag Models".
Lastly, for those not familiar with Mach numbers, I figured I'll do a formal introduction. When we say a bullet is travelling at Mach 2, we mean that its speed is twice the speed of sound, and so on and so forth. Writing supersonic speeds as Mach numbers makes things much cleaner-looking.
Atmospheric Density Variations
A sharp eye might notice that in our equations of motion so far, we'd taken air density to be a constant—and this might not be reasonable for rifle bullets that might travel kilometers into the sky. We expect air density to diminish with altitude, and perhaps we ought to account for this properly.
Here I reference an examination question as given in the 39th International Physics Olympiad, Vietnam, Theoretical Problem 3. Titled "Change of Air Temperature with Altitude, Atmospheric Stability and Air Pollution", the problem provides great insights into the functioning of our atmosphere. The first few sub-problems on temperature and pressure variations fit our needs perfectly.
As long as we stay within the troposphere, it is reasonable to model atmospheric temperature as a linearly decreasing function of altitude. The troposphere is 6-10 kilometers high, and we fully expect rifle projectiles to stay within this layer of the atmosphere; and so this is perfectly fine for our purposes.
Above we have introduced Λ, the temperature lapse rate of the atmosphere. Its empirically characterized value is about 0.0065 K/m, as stated for the International Standard Atmosphere. Assuming that the atmosphere is in a state of hydrostatic equilibrium, and that air is an ideal gas, we can solve for atmospheric pressure. Consistent with the usual conventions, M represents molar mass and R represents the ideal gas constant.
And from here it is straightforward to obtain atmospheric air density.
We can plot our calculated atmospheric air density out, and we can see that at a kilometer above ground—a height easily attainable by rifle bullets—the air density already drops to 90% of its sea level value.
The modified equations of motion with quadratic drag would then look something like this.
And here comes the million-dollar question: does it matter? And the answer is that it does, especially at higher launch angles. If we plug in reasonable values for the various parameters—mass on the order of a few grams, drag coefficient of around 0.30, muzzle velocity in excess of 900 meters per second, standard atmospheric temperature and density—we find that there's a considerable difference between our model accounting for air density variations, and the original that assumes constant air density.
And so air density variation is indeed a significant factor. We'll be using the modified equations of motions for analysis henceforth.
Bullet and Drag Models
Wikipedia lists the muzzle velocity of the SAR 21 rifle, the main workhorse of the military here in my home country, to be 970 m/s with the M193 cartridge—which is probably, and this is largely a guess on my part, in use in modern times. And so from here on we shall take the muzzle velocity to be 970 m/s, and the cartridge in consideration shall be the M193. It looks something like this, in comparison with the older M855 cartridge.
And here we have to, of course, remember that the bullet is actually a really tiny portion of the entire cartridge—it's the cone-shaped piece at the very tip. What comes out of the barrel is the bullet, and only the bullet; and so everything that we are going to consider—the mass of the projectile, the shape, the aerodynamic properties—are all of the bullet itself. The remainder of the cartridge is irrelevant.
What we really need is a detailed schematic of the M193 bullet. Fortunately a factory drawing of the M193 is available online, and this allows us to model its geometry with relative ease. Do note that measurements in the drawing are in inches.
In addition, we also learn that the projectile weight of the M193 cartridge is about 55 grains, which translates to about 3.56 grams. We shall take its aerodynamic cross-sectional area to be its frontal area; its radius is easily gleaned from the schematic.
McCoy and Kolbe Ballistics Calculator
The first drag model we'll consider is one that we have reasonable confidence in—a program written by a researcher in the field of external ballistics, which had some forty years of history and numerous revisions. Documentation has it that it was first written in Fortran, and now a Perl version of it had been coupled to a graphical front-end and put up on the web, which we will be using. The expected standard errors in its predictions—as benchmarked against some rather extensive experimentation—are stated to be 3% in the supersonic region, 11% in the transonic region and about 6% in the subsonic region. We'll denote this as the Kolbe model, after the author of the Perl revision.
Entering the geometry of the M193 projectile into the program, we obtain the following graph of predicted drag coefficient.
Trouble is, the program does not provide a tabulated version of these drag coefficients, which we need for our analysis. So we employ some Mathematica skills to extract our desired information. We select only the pixels that are red in colour, transform their coordinates to reflect the axes scales, and then fit a Bezier function through them to average out discretization errors. I'll leave the technical details out of this post.
We end up with an extracted drag coefficient function, which looks something like this when plotted. The sharp peaks have been rounded slightly, but the error introduced by that would likely be negligible.
So this is the Kolbe model. What else is there?
In the worlds of firearms and ammunition, it is in reality much more common to refer to the aerodynamic characteristics of bullets using the ballistics coefficient, rather than the drag coefficient. The ballistics coefficient is very similar to the drag coefficient, but it is defined with respect to a reference projectile—it is the ratio of the projectile's drag characteristics against that of the reference.
The reason for doing this is convenience, at least in manufacturing and in commercial trade. We'd seen the non-linearity of the drag coefficient, especially in transonic regions; people figured we can bury this non-linearity by defining a coefficient with respect to a reference projectile, since the general shape of the drag profile is roughly the same for all projectiles. This allows them to communicate a single number, instead of an entire graph. And there are numerous reference projectiles defined to date, to accommodate the shapes of different projectiles—the most common ones used nowadays are the G1 and G7 standard projectiles, with the latter much closer in geometry to the M193.
Accompanying each of these reference projectiles is a reference drag coefficient profile. Ammunition manufacturer sites offer these reference charts; we obtain ours from Barnes Bullets LLC.
We can employ the same programmatic image-analysis technique and extract the ballistics coefficient profiles out of these images. And there's also a website offering G1 ballistic coefficients mined from commercial sources, and G7 coefficients mined from military sources—we learn that the G1 coefficient of the M193 projectile is 0.227, and the G7 coefficient is 0.122. So this allows us to convert everything to drag coefficient.
We'll denote these two drag models Barnes G1 and Barnes G7.
It is quite a pity, because I can't seem to find papers with experimentally-measured drag coefficients for the M193 projectile. I did find some published data on the M855 and SS-109 projectiles, though. They're in a report titled "Aerodynamic and flight dynamic characteristics of the new family of 5.56 mm NATO ammunition", published by the US Army Materiel Command.
We won't be using these drag profiles in our analysis, because they're for different cartridges altogether; but they are good examples to illustrate a point brought up earlier. It was explained that people use ballistic coefficients in place of the drag coefficient, in order to reduce an entire graph to a single number—and we can indeed observe that the drag profiles for different cartridges are of the same rough shape. This similarity is what allows this method to be feasible.
If we were to look at the profile of a 7.62 mm projectile, or taken to the extreme, a 0.50 caliber one, we would find that this similarity diminishes. In fact, if we were to change the geometry of the bullet, say, take away the boat-tail, or modify its spitzer curvature, we would ruin the similarity. This is the reason for having numerous reference projectiles, from G1 to G7 and beyond—so that at least one of these references would have a similar geometry to the projectile studied, and the error incurred in using the ballistics coefficient representation would be acceptably low.
Computational Fluid Dynamics Simulation
So far we'd seen numerous drag coefficient profiles, as obtained from literature and ballistics calculators. A questions arises—can we compute the drag coefficient profile of the M193 projectile ourselves, fully accounting for the geometry of the bullet and the physics of supersonic flight?
And indeed we can. This discussion would span multiple subsections, but I guarantee it's going to be exciting and fun and worthwhile.
Fluid Dynamics and Heat Transfer Physics
Let's take a look at the physics first. The primary governing equation for fluid flow is the Navier-Stokes equation, which can be derived from energy conservation and Newton's Laws, as applied to fluids. We also have the continuity equation, which establishes conservation of mass—in other words, making sure that portions of fluid doesn't just disappear into thin air. And we'll also need to account for the thermodynamics of the fluid flow, especially since we're potentially in the supersonic regime; this necessities a heat balance equation. And along with these differential equations, we also have a no-slip boundary condition on the bullet surface.
The derivation of the Navier-Stokes equation and the heat balance equation is complicated, and we will not examine them here—we have, however, written them in terms of the viscous stress tensor τ, because we'll be referencing it later on in our analysis.
Above u represents the fluid velocity, ρ represents the fluid density, p is the fluid pressure, T is the fluid temperature, S is the strain-rate tensor, ∇ is the del operator, and colons represent the double inner product. They might look complicated, but they're a set of standard equations that modern simulation softwares solve for, and what I have done is to merely drag them out from documentation.
What do the equations physically signify? We can dissect the Navier-Stokes equation first. In brief, the left-hand side represents inertial forces; and the first term on the right-hand side represents pressure forces, and the second term, viscous forces, and the third term, external forces. This single equation unites these various factors in the description of fluid flow.
The continuity equation is relatively straightforward. It says the divergence in the flow is balanced by changes in density, which makes intuitive sense. And the heat balance equation says that changes in fluid internal energy at a point is caused by thermal transfer—the different possible types of it all summarized in a single heat flux vector q—and viscous dissipation, and work done via compression, and external heat sources Q.
This is not all, I think—and here we're going into territory that I have a shallow understanding in. I had come across this whilst writing up this post, and I thought it might be worthwhile to share. The continuity equation seems to be a hyperbolic partial differential equation; and the Navier-Stokes and heat balance equations are elliptic for stationary flow. If we assume diffusive effects to be negligible, which is typically acceptable in high-speed flows, all three equations become hyperbolic. And a hyperbolic equation has the property that disturbances propagate outwards at finite speed—which is the emergent mathematical analogue of the finite speed of sound in fluids. Neat!
Of course, we cannot expect the thermal conductivity and viscosity of air to be a constant across the wide range of temperatures involved. This necessitates two additional equations, stemming from Sutherland's law.
Sutherland's law is derived from the kinetic theory of gases, with an idealized intermolecular potential. It provides an excellent description of the changes of thermal conductivity k and viscosity μ of the fluid, with respect to a reference value and a reference temperature. A fluid-specific constant, termed Sutherland constant, is also introduced.
We'll concern ourselves with the physics simulation of a large, but finite, volume of air around the bullet. For obvious reasons we would not be able to simulate an infinite air volume, or anything close in size to the actual atmosphere; we'll elaborate on the simulation method and computational limitations in the next subsection. Simulating this finite volume of air also means that we'll need boundary conditions for remaining boundaries of the volume, which we present below.
The idea is that instead of simulating a moving bullet in stationary air, we'll simulate a stationary bullet in moving air—these two should be equivalent. Such an approach is typically known as a virtual wind tunnel, and had indeed seen use in external ballistics. And this necessitates an air inlet and outlet, which is to be placed a significant distance in front and behind the bullet, respectively. We'll be exploiting axial symmetry—this is elaborated upon in the next subsection—and we therefore also have a single circumferential boundary.
Of all the boundary conditions above, most are pretty self-explanatory, except possibly the last one. It is a no-penetration condition, which dictates that no streamlines are to cross the circumferential boundary; and in addition we also dictate that there are no viscous effects at the boundary, which must of course be the case, for the boundary doesn't physically exist.
When faced with such a system of partial differential equations and complex geometry, analytical solutions are almost always unobtainable. We have to turn to numerical solutions; and one way of computing them is through the finite-element method, which had seen extensive use in simulations across numerous areas of physics, including computational fluid dynamics.
The finite-element method relies on geometry discretization—the splitting of geometry domains into a tessellation of small cells, typically triangles and quadrilaterals for 2D geometry, and tetrahedrons and cuboids for 3D ones. The collective tessellation is known as a mesh. The physics is then solved across the individual elements of the mesh; and because the elements are small and are of simple geometry, it becomes possible to make certain simplifications—such as linearization of the equations, or presuming that basis functions are polynomials. Mathematically, this process involves linear algebra with matrices of epic sizes.
The idea is that the numerical solution would improve in accuracy as the discretization is made finer and finer; and therefore a solution of any quality can be obtained, so long as one has sufficient computational resources.
Here comes our first problem. Our simulation has to cover both subsonic and supersonic regimes; and the thickness of a shock in supersonic flows is in the sub-micron range, which is really small as compared to the scale of the bullet. We can exploit the axial symmetry of our system, which constrains us to a zero-yaw assumption; but even so, we would require at least several hundred million mesh elements to capture such fine details. It is clearly not feasible to solve for several hundred million degrees of freedom.
And so, in addition to exploiting symmetry, we also utilize adaptive mesh refinement. Instead of using a very fine mesh from the start, we start with a coarse one, and iteratively refine it where needed. More specifically, we refine the mesh where the error norms are large, which typically corresponds to regions where there are changes in the solution fields too sharp to be captured by the mesh. This is, in fact, an approach taken by SpaceX in their numerical re-entry and launch simulations.
We can see the mesh refinement in effect in the following mesh plots. The top row was extracted from a subsonic flight computation at Mach 0.5, and the bottom row was extracted from a supersonic computation, at Mach 2. From left to right, we present the initial mesh, first refinement, second refinement, and third refinement. We can see that for subsonic flow, the refinement occurs mostly in proximity to the bullet, where the boundary layer is; and for supersonic flow the refinement tends to follow the shocks.
With this approach, we can capture the necessary details in the fluid flow whilst maintaining an overall low mesh size; it thus becomes feasible to run the simulation on a personal desktop. Numerical convergence remains a huge issue, especially in supersonic regimes—and it is rather inevitable that extensive manual tweaking has to be involved. I must explain here that this tweaking does not refer to modifying the underlying physics in any sense; what is tweaked are the numerical solvers and their algorithmic parameters.
Another thing—we would take a leap of faith and assume that the fluid flow is laminar. Yes, it is true that the Reynolds number suggests very strongly that there is turbulence in our velocity regimes, and yes, literature does indicate that turbulence is certainly present in trailing regions. And so we should ideally be accounting for turbulence. Throwing turbulence into the mix really does mess with the convergence, however, and I had been trying the k-ε and k–ω turbulence models, to no avail. I ended up settling with laminar flow, and besides—do take this with a pinch of salt—I think there wouldn't be a significant difference in drag force results. In my defense our simulation results do seem to match closely with the experimentally-verified Kolbe model, as we'll see later, so maybe this simplification isn't too costly after all.
A selection of simulation results are presented below, covering both subsonic and supersonic flight. A more extensive set of results can be downloaded at this link, for those interested.
—— Fluid Mach Number Plots ——
—— Fluid Pressure Plots ——
—— Fluid Temperature Plots ——
Let's have a look at the top left simulation, corresponding to a subsonic projectile velocity of Mach 0.45, first. We observe a localized high-pressure, low-velocity spot at the tip of the bullet, which is fully expected as the bullet rams into the stationary air ahead. Near the boat-tail we observe flow separation, and immediately behind the bullet is a region of significant rarefaction—this is where we would expect turbulent effects, such as vortex shedding, to be most prominent, if we had accounted for turbulence in the simulation.
Very interestingly, we can observe with great clarity the formation of shock waves as the bullet crosses the sound barrier. We first observe the formation of normal shocks, as the deflected airflow around the bullet exceeds the local speed of sound—this is very similar to the normal shocks formed on the wings of a transonic aircraft. At slightly higher velocities, these normal shocks can no longer keep up with the bullet, and they coalesce into a large trailing shock. And at even higher supersonic velocities, a powerful bow shock forms, with near-discontinuous changes in air velocity, pressure, and temperature.
So far we had been discussing shocks from a taken-for-granted perspective, without actually touching upon why they form in the first place. They arise because the propagation of disturbances in any fluid is fundamentally limited by the speed of sound in that medium—and this must be true, because sound is nothing but vibrational motion passed from one molecule to another, so by definition transfer of disturbances cannot happen faster than its speed.
We can think of it this way: when air flow exceeds the speed of sound, air molecules in that region would have no way of knowing what lies ahead, for information—the changes in pressure that tells the molecules how to act—cannot reach them fast enough. And they therefore pile up somewhere, like a traffic jam, and they do so with sufficient extent to create a localized, almost-discontinuous jump in gas properties, which is what we term a shock.
Digression #3: Mach Angle
This will be a short intermission, seeing that we're already talking extensively about shocks . We had said that shocks form because of the finite speed of sound; and we should therefore expect a rather simple geometrical relation between the shock angle, and the speed of sound and speed of projectile. This is well-known in physics, with the shock angle termed the Mach angle, and the axially symmetric shock wavefront termed the Mach cone.
We can attempt to verify this relation using our simulation results. We can sample the pressure changes as given by our simulation, say, on a 10000 by 10000 grid surrounding our bullet. We then apply a gradient filter to this scalar grid, filter for the outermost portions of the shock, and perform a linear regression to capture its angle. All these can be done computationally without much trouble in Mathematica.
And we can indeed see that our simulation results matches our expected analytical Mach angle expression perfectly.
The reason why we sample only the outer regions of the shock cone, is that near the bullet the shock curves—and this is indeed the ambiguity we had mentioned earlier in "Reynolds Number and Quadratic Drag". Near the bullet the temperature and density of the air changes so much that the speed of sound is significantly modified; and this in turn means that the Mach angle would change, producing the curvature seen. Because we had taken sea-level speed of sound as a reference, it is appropriate for us to sample the shock far away from the bullet, where the gas properties are closest to atmospheric.
Simulation Drag Coefficient Results
Let's get back to what we're interested in—the drag coefficient of the bullet. We can easily re-arrange the drag equation for this purpose, expressing it in terms of the drag force.
In general, drag force can be decomposed into two components—viscous drag, and a pressure-induced force. The former arises because portions of the fluid close to the bullet tends to get dragged along by the bullet surface, and the viscosity of the fluid then transmits this shear force to nearby fluid layers; the result being that the bullet experiences a resistive viscous force. The pressure-induced drag, separately, arises because the movement of the bullet through air tends to compress the air in front of it, and expand the air behind it. This pressure gradient tends to push the bullet backwards, resulting in a drag force.
We can evaluate the viscous drag force by integrating the axial component of the viscous stress tensor across the bullet surface. Similarly, we can evaluate the pressure drag force by integrating the air pressure, projected onto the bullet surface.
And here we present our final drag coefficient results. To obtain the best-fit curve, the raw results were first smoothed through a moving average, and then spliced using a Bezier function. We observe a sharp peak in the transonic region, similar to the other drag profiles; and we also observe a gradual decrease in drag coefficient in the supersonic regime.
In "Reynolds Number and Quadratic Drag", we had discussed that the relative importance of viscous effects diminishes as projectile speed increases. And we can indeed see that this is the case, from our simulation results. As the projectile speed increases, the percentage contribution towards total drag force of pressure-induced drag increases steadily, whereas that of viscous drag decreases. We see a slight jump as the bullet crosses the sound barrier, and in the supersonic regime viscous drag is almost entirely negligible.
Projectile Lethality and Terminal Ballistics
Summary of Drag Models
Here's where we are going to finally answer our question. Let us first summarize all the drag profiles we're considering in a single graph. Amongst them are the Kolbe model, the Barnes G1 and Barnes G7 models, and the drag profile we had computed via finite-element computational fluid dynamics (CFD) simulation earlier. We can observe that the CFD and Kolbe profiles are very similar, differing slightly only at the transonic peak and the supersonic tapering—this is a good indication that the simulation is reliable.
The constant coefficient model is derived from the velocity-averaged Kolbe profile, as a benchmark baseline. We would use this to see if a constant coefficient approximation is reasonable, after all.
We can plug these drag profiles into the equations of motion we had derived earlier, and plot out the ballistic trajectories. Below we present a graph of these trajectories, for arbitrarily-chosen firing angles of 35, 70 and 85 degrees. As expected, we can observe that the CFD and Kolbe models yield very similar results, and that the Barnes G1 and Barnes G7 are mutually similar.
To gauge the lethality of the projectiles, we would have to understand the behaviour of the projectiles as they impact the human body. Such a field is known as terminal ballistics, and we would dedicate the next subsection to a brief discussion of what's relevant to us.
The lethality of the projectiles depend on a lot of factors; and one of them is the bullet type. For instance, hollow point projectiles are designed to allow soft matter to flow into its frontal hollow, and in doing so, expand its diameter, so that it disrupts more tissue; this has the effect of creating more serious wounds at the cost of penetration ability. Full-metal jacket ammunition, on the other hand, does not have have this expansion effect, and creates narrower wound channels with greater penetration ability. There are many other designs, such as ballistic tips and soft points, which attempt to achieve some sort of enhanced middle ground between hollow points and metal jackets.
Typically ball cartridges are used for standard-issue military rifles; and the M193 ball cartridge is a metal-jacketed, lead alloy core bullet. It is interesting to note that hollow-point projectiles are strictly prohibited for military purposes under the Hague Convention of 1899 as they cause excessive wounding. And so our discussion shall focus on the metal-jacketed M193 ball projectile.
Another major factor, of course, is the impact velocity of the projectile. We should understand that impact velocity affects not only the depth of the wound, but also the type of tissue damage.
I am neither a biologist nor a terminal ballistics expert; and so I shall attempt to offer only a shallow, qualitative overview. The bulk of the information presented here was sourced from published literature, such as "Gunshot wounds: A review of ballistics related to penetrating trauma" by Stefanopoulos et al., "An Analysis of the Wounding Factors of Four Different Shapes of Fragments" by Ma et al., and "Comparison of the Terminal Ballistics of Full Metal Jacket 7.62-mm M80 (NATO) and 5.56-mm M193 Military Bullets: A Study in Ordnance Gelatin" by Ragsdale et al. They contain detailed investigations that would explain terminal ballistics far better than I can.
For low-velocity subsonic impacts, the primary wounding mechanism is tissue crushing—that is, the overpressure directly in front of the penetrating projectile squeezes and crushes bodily tissue. The characteristic for such a type of injury is a narrow wound channel, barely larger than the projectile diameter; and because of the relatively small size of the wounded region, such an injury is unlikely to be fatal, unless it strikes the central nervous system and causes nervous shutdown, or major arteries with consequent lethal blood loss.
At velocities of 600 m/s or greater, cavitation becomes important. The high-velocity projectile generates an expanding bubble in trailing regions as it moves through soft tissue; and this bubble can easily expand to numerous times the diameter of the projectile, shearing and tearing tissue, even breaking bones indirectly. The damage caused by cavitation can be extensive, and it increases the lethality of the projectile significantly.
Perhaps the most lethal effect of all is that of fragmentation, where the impact is vigorous enough to break the projectile into numerous separate shards, each of which creating its own wound channel and cavitation bubble. Most modern projectiles are designed to achieve considerable fragmentation, to maximize effectiveness. Fragmentation requires a minimum impact velocity, and also a minimum distance to occur; it seems the threshold velocity for the M193 ball projectile is about 800 m/s. There is literature pointing towards delayed fragmentation of the M193 in certain cases—that it begins to fragment only after it exits the body—but we shall assume that effective fragmentation is achieved.
And so we can split the velocity range into four domains. In the subsonic regime, lethality is unlikely; and in the Mach 1 to 600 m/s regime, wounds are severe, and lethality is possible. From 600 m/s to 800 m/s lethality is likely due to cavitation; and beyond 800 m/s lethality is highly likely, due to fragmentation. We shall refer to these as the subsonic, intermediate-velocity, high-velocity, and high-supersonic regimes respectively, consistent with literature.
It should also be mentioned that the survivability of a gunshot wound is also highly dependent on the quality and promptness of medical aid administered. We have essentially assumed a level of medical intervention identical to that considered in literature. More advanced medical aid could shift the lethal impact thresholds higher than what we have considered here.
Projectile Impact Velocity
Now we evaluate the expected impact velocities of our projectiles. Let us assume that the projectile impacts occur at ground-level, as it must be, if it is to strike a pedestrian on the streets. The system of equations for the impact velocity can be written as follows.
And we observe that for most firing angles, the M193 projectile would be travelling far too slowly by the time it comes back down to cause any considerable damage. In fact their terminal velocities seem to hang around the 130 m/s region, which is far less than the speed of sound. We can therefore conclude that for the case of ground-level impact, it is overwhelmingly likely that the bullet would be non-lethal.
It would be interesting to see in closer detail what really happens at low firing angles. It does seem that below a certain firing angle, the bullet would retain sufficient velocity to cause significant damage.
We can observe that the bullet enters the subsonic regime at angles above 0.5 degrees or so; and a high-velocity impact, which would likely be lethal, would require a firing angle of less than 0.3 degrees. This means one would have to aim pretty much directly at an innocent passer-by, in order to have a fatality incident via negligent firing—and I really doubt anybody would be as careless as that.
This, of course, applies only to ground-level impacts. It is perfectly possible that the unsuspecting victim happens not to be at ground-level—he could be in a building and the projectile happens to pass through a window and strike him, or he could be standing on a bridge, or he could simply be on elevated ground. Not restricting our impact calculations to ground-level makes things a little more interesting.
For starters, we now have two cases to consider: the first case being an impact while the projectile is still ascending, and the other while the projectile is descending. We therefore modify our equations, with h being the arbitrary height of impact.
We can evaluate the impact velocities at various heights; and below we present a set of graphs to illustrate the results. Click on any of them to see it in detail.
And the conclusion we come to for above-ground impacts is that regardless of the impact height and firing angle, the projectile's velocity during descend is almost always too low for lethality. Things are a little different for the ascending branch—we see that at impact heights of 100 meters or below, which is typical of modern residential buildings and bridges, the majority of firing angles do yield high-velocity projectiles capable of considerable harm. At higher elevations it becomes increasingly difficult to achieve lethal impact velocities on the ascending branch, and above 800 meters or so lethality is impossible.
And so we have, at last, answered our question.
The bullet can kill, but only if it is aimed almost directly at a target, or if it strikes a low-lying target during ascend. In all other cases lethality might be rather far-fetched.
There is, of course, another way to visualize the impact velocity results—we can instead plot the full velocity profiles of fired projectiles, against height attained. We mark the ascending branch as solid lines, and the descending branch as dashed lines. This allows us to see at what height exactly do the projectiles fall past the various velocity thresholds.
Of course we'll arrive at the some conclusion; this is just a different way of seeing things. The descending branch is almost never lethal, except at very low firing angles of a degree or less. On the ascending branch, the projectile can be lethal, so long as it is fired at a sufficiently high angle. The higher the target, the higher this angle needs to be.
The velocity profiles also gives us a neat way of visualizing how the projectile behaves along its trajectory. For instance, we see that at low firing angles, the projectile velocity is a smooth curve—the combined effect of air drag and gravity is to decelerate the projectile consistently. However, as the firing angle increases, the bend in the velocity profile gets sharper and sharper; this is consequent of the opposing effects of gravity and air drag on the descending branch, the former attempting to accelerate the projectile whilst the latter attempting to decelerate.
At very high firing angles, the projectiles slows down near the peak of its trajectory, and speeds up as it descends—this appears as a kink on the graphs. And at a firing angle of 90 degrees, the projectile comes to a complete standstill at its peak.
Probability of Lethal Impacts
We now know that lethal impacts are not very likely, unless the firer was aiming directly at a target—which shouldn't happen by mistake. This prompts us to ask a follow-up question: can we calculate the rough probability of lethally hitting someone, if one was to fire by mistake?
First of all, a disclaimer—we are concerned only with a rough probability, and by rough, we mean something similar to an order-of-magnitude estimate. The calculations presented in the following subsections are undoubtedly the least reliable and most dubious of everything written in this article; but even with such potential for inaccuracy, I still feel they are worth the space on this post, because tied together with those calculations is a surprisingly large amount of interesting sideline ventures. We'll be visiting, amongst others, computational geometry, some data mining, geospatial calculations, and a fair amount of mathematical modelling and number crunching.
There will be many approximation symbols along the way, and our final results are to be taken with a pinch of salt. But I hope this end-game adventure would be as interesting to you, as it was to me.
Heuristics and Mathematics
Quite obviously calculating the probability of lethally hitting someone is not a trivial task. There are too many factors involved. And so we establish some heuristics to help us with building the mathematics:
Impacts can be split into two primary cases—an above-ground impact in buildings, and a near ground-level impact on the streets and remaining land area.
Lethal impacts in buildings can only happen if the projectile enters through a window or glass facade. Structural concrete is assumed to present sufficient impedance to render projectiles non-lethal after penetration.
A projectile can at most impact once. A projectile cannot pass through a target completely and strike another lethally. The chance of a lethal impact after repeated penetration of windows or other building structural weak-points is negligible.
We shall denote the total probability of lethal impact as P, which would, of course, be a function of the firing angle. The first heuristic gives us the following decomposition of P into substituent probabilities, which we will treat separately.
Now we posit that the probability of lethal impact near ground-level can be given by the following expression. Amongst other factors, we have included in our consideration the areal density of the human population in our city, the areal density of streets and non-buildings, and the impact velocity of the projectile. A more in-depth explanation would be given shortly.
It is best to analyze our proposed equation term-by-term, so that we know exactly what it means. The first term is the areal density of streets and non-buildings in Singapore, and it is equivalent to the probability of the projectile landing on land area not occupied by buildings. This is a necessary condition for a near ground-level landing. Do note that what we're really referring to here—and in many other instances within the next few subsections—is the average areal density as evaluated across the whole of Singapore. Of course the areal density of buildings would vary from neighborhood to neighborhood and from district to district, but we can simplify things by taking the average, courtesy of our homogeneity and isotropicity heuristic.
The second term, as stated in the equations, is the projected areal density of the ground-level human population present here in our city. We say projected, because the areal density of our human population exposed to the projectile would vary with its incoming angle; and we say ground-level population, because at any point in time some percentage of our population would be indoors, shielded from ground-level descending impacts. We will examine this term in closer detail in a subsection below, "Projected Population Areal Density".
The third term is coined the lethal impact velocity threshold function for good reason. In our analysis earlier we had said that lethality is only likely when the projectile is travelling at supersonic speeds. Specifically, we had concluded that lethality was possible with impact velocities between Mach 1 and 600 m/s; anything above that, and lethality was very likely. We would simplify this and say that lethality is certain at supersonic speeds, and negligible at subsonic speeds. Mathematically this is equivalent to δ being a piecewise-constant function.
And the fourth term is the time-averaged probability that the projectile hadn't impacted anything up till its landing. This is to enforce our single-impact heuristic. We will explain this term in greater detail slightly later, because it requires some additional information to construct.
Moving on, the probability of lethal impacts in buildings can be written in a similar form, as presented below. It is necessary to introduce time into the equations, so that the probability of lethal impact at any instance during the projectile's flight can be calculated; we will get rid of this time dependence through an averaging process later on.
Like before, we'll break it down term-by-term to understand what it means. Similar to the expression for the near ground-level lethal impact probability, the first term is the areal density of building footprints in our city, or equivalently, the probability of the projectile being in a building footprint. And the second term is the areal density of glass on the facades of buildings, or equivalently, the probability of the projectile entering through a glass surface. As explained before, these are average areal densities, as evaluated across Singapore.
The third term is the projected above-ground population areal density, and the fourth is the impact velocity threshold function. We had visited these before. The fifth term is new—is it the probability that the projectile is at an altitude low enough to intercept buildings. In other words, it is the complementary cumulative distribution function of the building heights here in Singapore. We will have to obtain the distribution of building heights somehow, and this is explained in a later subsection, "Building Height Distribution".
And the last term is the time-averaged probability that the projectile hasn't yet impacted anything; again, a manifestation of our single-impact heuristic. Here we finally see the from of Q—essentially the probability of the projectile intersecting a building along its trajectory.
Before we attempt to solve for the lethal impact probabilities, we'll need to figure out the values for the various constants and the form of the various functions. This is what we will do in the next few subsections.
Building Footprint Areal Density
We shall start with the areal density of buildings. The simplest way of obtaining such an estimate is to attempt to look up construction plans and city development reports, and there is indeed one, the 2010-2030 land use plan publication by the Ministry of National development, Singapore. The publication lays out the percentage land area allocated for housing, industry and commerce, parks and recreational facilities, amongst many other categories; but it does not tell us how much land is actually used for buildings in each category. And so it is largely useless, unfortunately.
The Urban Redevelopment Authority also has a comprehensive masterplan, but it suffers the same caveat. We cannot obtain accurate area figures from it.
My various attempts to find literature on building footprint statistics have failed miserably, and so I will resort to my own methods. The idea is that I'll download high-resolution local map data and analyze them computationally—and calculate for myself the total footprint of all buildings in Singapore. Obtaining areal density from that would then be trivial arithmetic.
The first step is to download the necessary geographical data. I had done so by grabbing one of the OpenStreetMap data packages from BBBike. OpenStreetMap is a crowd-sourced map of the world to which anyone can contribute; and being an open project, its data repositories are publicly available. Visualization and processing of the data was then done in QGIS, an open-source geographical information software.
The raw geographical data of Singapore, along with some neighbouring parts of Malaysia, looks something like this when plotted.
We can command the software to filter only for buildings. This gives us a much cleaner plot.
And because the building footprints are all polygons, calculating their area programatically is easy. We obtain a total building footprint of 84,178,390 square meters. Separately, the Department of Statistics Singapore lists our total land area to be 719.7 square kilometers. And so we can evaluate the footprint areal density of buildings rather easily.
For those interested in the details on how the geographical data was processed, and the calculations done, please read my follow-up blog post titled "Manipulating Geographical Data with OpenStreetMap and QGIS". The technical details are elaborated there, and other aspects of geospatial computation are also visited.
Building Height Distribution
Unfortunately OpenStreetMap does not offer data on building heights, and so we cannot employ the same approach as we had for building footprints. And, just as we couldn't find any literature on building footprints, we couldn't find any on height distribution as well.
We're in luck, because there's another avenue—Emporis, a real-estate data mining company with a huge database on infrastructures worldwide. They have a website available to the public, on which we can see limited data for about 20 buildings per page. There are more than 8000 buildings in Singapore. To obtain complete datasets conveniently one would have to make a purchase—but fear not, because here comes Mathematica.
We can easily write a program to visit all of Emporis's webpages and extract the building height data that we require. In fact, we can go one step further, and extract all the building data that the website has to offer, including buildings names, number of floors, building type, year of construction, amongst others. Some of these additional information would come in handy, as we'll see later.
With this, we are able to obtain quite an extensive measure on building heights here in our city. We can visualize the data by plotting histograms, one for the combined dataset, and one each for the various districts in Singapore.
We can combine these to create an overall probability density function, which we require for our calculations. In doing so we have converted a discrete dataset into a continuous distribution, and this involves some mathematics that we would not elaborate here. For those interested in the technical details on how the data extraction and processing was done, please read my follow-up blog post titled "Extracting Data from Websites with Mathematica".
Glass Areal Density on Building Facades
Unfortunately, I was not able to find concrete data on the window area on buildings, as I had done for building footprints and building heights. I had been pretty certain that there would be building codes or architectural guidelines on the total window area for high-rise buildings; but that intuition turned out wrong, for I couldn't find anything of semblance to a hard rule. Architects gauge suitable window placements and window areas by calculating the expected illuminance inside their buildings, and this is highly dependent on building geometry.
I have, however, read that a general rule-of-thumb is that window area should not exceed 25% of the total floor area of a house. If we apply this approximate guideline to the data that we at our disposal, we can indeed calculate the areal density of glass surfaces. Below we let Pr represent the perimeter of a polygonal building footprint, and h represent the height of a building. An overhead bar represents an arithmetic mean.
The mean building footprint area and the mean building perimeter can be obtained from our OpenStreetMap dataset; and the mean number of floors and mean building height can be obtained from our Emporis dataset. In this way the areal density of glass surfaces can be calculated.
Projected Population Areal Density
We had said earlier that the projected population areal density depends on the incoming angle of the projectile. This must, of course, be the case. Say, if the projectile was to approach from directly above, it can only impact the head and shoulders of a pedestrian, which is a rather small presented area; whereas, on the other hand, if it was to approach horizontally, the entire frontal area of a pedestrian would be exposed. The area of a human target differs, as seen from the bullet's perspective.
In physics we oftentimes simplify the geometry of complex objects to more manageable shapes. But why simplify anything, when we have Mathematica?
We can manipulate 3D models of the average human in Mathematica with ease. Calculating the exposed area of a human at varying projectile trajectory angles is equivalent to calculating a 2D projection of the 3D human body; and this can be done computationally. Below we present some rendered graphics of the human body as viewed from the perspective of an incoming projectile, and along with those, the projected 'shadow', at various elevation angles α and two different azimuthal angles β.
Like the other technically-intensive tasks, I would not go through the exact programmatic details and mathematics here. I have written a separate blog post about this, titled "The Human Body: 3D Image Projection in Mathematica", for those interested in further elaboration. It suffices to understand at this point that the projection and the area calculations can be done, and that high-resolution results can be obtained. A sample of these results are presented below.
Notice that we need only perform the computations over an 180-degree range of α and β, due to the intrinsic symmetry of the human body. We now recall our homogeneous and isotropic assumption—we can reason that we are applying our calculations to such a large land area that the individual orientations of potential human targets can be taken to be completely random and uniformly distributed. In that sense, we can define the average projected area of a human target to be the following.