Anyone who brought it upon themselves to write an overly realistic bullet simulation for their game will inevitably come to a point where they discover – this isn’t so realistic, why do bullets just disappear when they hit a wall?. While the majority of game developers shrug, research the subject for about 13.2 seconds, ignore the subject completely and wear a blank expression for the rest of their programming days, I, dear reader, decided to embrace the pain.
Panoptesv has published an excellent write-up on different penetration regimes [Accessed 9th of August, 2019]. This guide details pretty much every situation you can dream of in projectile penetration. Both non-deforming- and deforming bullets. Low, high and intermediate speed regimes – although these are only analytically solvable for the non-deforming bullets. There’s even some pointers on how to handle thin wall penetrations. It’s definitely not complete, but it’s a good beginners’ starting point compared to much of the published literature from the 70’s on the subject.
What is the problem?
The problem with Panoptesv‘s excellent article is that what is given is the penetration depth of these situations, but not this depth with time nor the remaining velocity. This is where the problem lies – the frame rate of your bullet simulation will likely not exceed about 60 FPS, or 16.6(7) ms per frame. This is plenty fast if your simulations occur completely within one medium – say, air. But, as soon as another medium, such as drywall, is entered, you are working with fractional physics simulation steps. An example:
The time frame of a 9x19mm Para bullet (diameter = 9 mm, mass = 8 g) penetrating something under the assumption of non-deforming bullets is not going to exceed about 1 ms, and it will penetrate into the target for about 88 cm. So? That means penetration is much faster than the framerate – why bother simulating it? Because, the residual velocity of the bullet depends on how deep your bullet penetrates. And, if your bullet penetrates 50 cm of drywall when it could penetrate 88 cm, it’s gonna come out the other end. You’ll need to know how far along the penetration your bullet is at that point, otherwise you can’t come up with accurate residual velocities. Furthermore, if your physics simulation has less than 1 ms left to be simulated, your bullet will need to stop being simulated in the middle of the drywall, only to be further simulated as the next physics simulation frame comes around.
How are you gonna solve it?
I’ve spent about a full day (13:00-21:00, a typical day) converting the non-deforming bullet formulae from penetration as a function of velocity to penetration as a function of time, and time as a function of penetration. I made several lucky guesses which sped up the process, but then I hit the roadblock that is called What is the average value of x(t) = Log(1+x'(t)^2)? I’ve solved it all, though, here’s the results (and a nice graph to show you I solved something! Fancy!).

Non-deforming projectiles: Low speed penetration
I’ll give it to you straight, like a pear cider made from 100% pears. Here’s the answers you’re looking for, the rest of this section (1 out of 3) is the derivation. There’s three sections in total, each built up much like this one. Feel free to skip the derivations and straight up copy the final results at the start of each section! I know I would.
Penetration Depth x(t) [m] as Function of Time t [s]
x(t) = t * v(0) - A * Yc * t^2 / (2 * m)
const = m / (2 * A * Yc)
x(t) = t * v(0) - t^2 / (4 * const)
Penetration Time t [s] as Function of Penetration Depth x [m]
t(x) = 2 * const * v(0) - sqrt(4 * const * (const * v(0)^2 - x))
or..
t(x) = 2 * x / (v(0) + sqrt(v(0)^2 - x / const))
Now let’s begin the derivation.
x = m * vi^2 / (2 A Yc)
x = const * vi^2
For this case, let us define a constant const = m / (2 * A * Yc) for brevity, with m [kg], the mass of the projectile, A [m^2], the frontal surface area which could be calculated as A = PI (d/2)^2 with d [m] the bullet diameter, Yc [Pa], the cavitation strength of the target which apparently matches about three times the yield strength of a material. Furthermore, we have x [m], the penetration depth at a certain time, and vi [m/s], the velocity upon impact.
x(vi) = const * vi^2
We should immediately notice that Panoptesv did us dirty – this entire formula is not time-dependent! How can we possibly know what the velocity is at a certain penetration depth, for example? Calm down! I randomly guessed the following, and miraculously it worked. How do you know it worked? So I tested the formula by checking that the distance x(vi = 500, v = 250) plus the distance x(250, 0) indeed matches x(500, 0). Later I checked that x(600, 599) + x(599, 598) + (..) + x(3, 2) + x(2, 1) matched x(600, 1). Again it did – this was a lucky guess.
x(vi, v) = const * (vi^2 - v^2)
So! That was easy! What’s next? Well, we only know the distance travelled between velocities vi and v. We’d like to also know the time it takes between those two velocities. My first instinct was to over complicate it, so instead I looked at the simplest way to do this. The distance travelled can also be described as x = avg(vi, v)*t with avg() a function of both velocities and possibly more, and t [s] the time in seconds from the start of the penetration. What average do we want to use, though? I guessed correctly that it was the arithmetic mean. First try.
x = avg(vi, v) * t
t = x / avg(vi, v)
(..guessing..)
t(vi, v) = 2 * x / (vi + v)
Again I tested my typical scheme – t(500, 250) + t(250, 0) did indeed match t(500, 0); great success! This is a problematic result, though, because this equation has three unknowns, both x and v are values we don’t know and their resulting value t is similarly unknown. What now?
Luckily, we also have our previous equation for x, which can be rewritten to read v = sqrt(vi^2 - x/const), which can be filled in to yield t(vi, x).
x = const * (vi^2 - v^2)
vi^2 - v^2 = x/const
v^2 = vi^2 - x/const
v = sqrt(vi^2 - x/const)
t = 2 * x / (vi + v)
t = 2 * x / (vi + sqrt(vi^2 - x/const))
While this is a nice way to calculate the time taken to reach distance x with starting velocity vi, it will be quite convenient as well to know the maximum distance x we can penetrate within our remaining time t. Let’s do it the only way we know how: calculus.
t = 2 * x / (vi + sqrt(vi^2 - x/const))
x = t (vi + sqrt(vi^2 - x/const)) / 2
x = t * vi / 2 + t / 2 * sqrt(vi^2 - x/const)
And then we’re kinda stuck. Luckily, WolframAlpha comes to our rescue:
x = t * (4 * const * vi - t) / (4 * const)
x(vi, t) = t * vi - t^2 / (4 * const)
Quaint!
Wait.. did I just say WolframAlpha? Maybe I made a mistake somewhere… Better check absolutely everything again but only using WolframAlpha.. right? The only way we can convert our (x, v)-based formulae to (x, t) formulae is through the fact that v = dx/dt = x'(t), e.g by forming a differential equation. Boooooo.
x(t) = const * (vi^2 - v^2)
v = dx/dt = x'
x(t)/const = vi^2 - x'(t)^2
x'(t) = sqrt(vi^2 - x(t)/const)
I haven’t a clue how to solve this. But. WolframAlpha does. And the result is.. Urk. Luckily this gross thing can be rewritten a bit more neatly by spotting -(2 c1 t + c1^2 + t^2) = -(c1 + t)^2.
x(t) = (- 2 * c1 * t - c1^2 + 4 * const^2 * vi^2 - t^2) / (4 * const)
x(t) = (4 * const^2 * vi^2 - (c1 + t)^2) / (4 * const)
x(t) = const * vi^2 - (c1 + t)^2 / (4 * const)
rewrite to t
t(vi, x) = c1 - sqrt(4 * const * (const * vi^2 - x))
So, I haven’t a clue what’s the value of c1 which was added as a result of solving the differential equation, somehow, I forgot how it all works. But. When I compare the results if t from WolframAlpha with the results I got without it, it seems the equation starts the wrong way around – when t(vi, 0) = tmax, while t(vi, xmax) = 0. Huh. I guessed that c1 = tmax, because this would reverse the order again.. fixing the problem. Filling in c1 = t(x=0), everything works.
c1 = tmax = t(vi, 0)
c1 = sqrt(4 * const * (const * vi^2))
c1 = sqrt(2^2 * const^2 * vi^2)
c1 = 2 * const * vi
filling in for t
t(vi, x) = 2 * const * vi - sqrt(4 * const * (const * vi^2 - x))
This result is somehow different from the other, but they give the same values. Finally, this whole ordeal can be rewritten to x(vi, t).
t = 2 * const * vi - sqrt(4 * const * (const * vi^2 - x))
(t - 2 * const * vi)^2 = 4 * const^2 * vi^2 - 4 * const * x
4 * const * x = 4 * const^2 * vi^2 - (t - 2 * const * vi)^2
4 * const * x = (4 const^2 vi^2) - (4 const^2 vi^2) - t^2 + 4 * const * vi
4 * const * x = 4 * const * vi * t - t^2
x(vi, t) = t * vi - t^2 / (4 * const)
Non-deforming projectiles: High speed penetration
Penetration Depth xp(t) [m] as Function of Time t [s] (v >> vthr)
xp(t) = 2 * m / (Cd rhot A) * ln(-(-2 * m / (Cd rhot A) - t * vi) / (2 * m / (Cd rhot A)))
xc = m / (Cd rhot A)
vthr = sqrt(2 * Yc / (Cd rhot))
xp(t) = 2 * xc * ln(-(-2 * xc - t * vi) / (2 * xc))
Penetration Time t(xp) [s] as Function of Depth xp [m] (v >> vthr)
t(xp) = 2 * xc / vi * (exp(xp / (2 * xc)) - 1)
Next up is the derivation of high speed penetration. Snicker.
xp = 2 * xc * ln(vi / vthr)
xc = m / (Cd rhot A)
vthr = sqrt(2 * Yc / (Cd rhot))
What is very different about this case is that the penetration depth xp listed is the depth that is reached until v = vthr, with vthr [m/s] a threshold value below which the simplification of high speed penetration is no longer valid. xc [m] is a length scale, calculated from some constants. Both depend on Cd [.], the drag coefficient of the bullet upon penetration (typically Cd = 1.0 in this situation, according to Panoptesv), rhot [kg/m^3], the density of the target (e.g rhot = 1000 kg/m^3 for water). All other variables have already been explained for the low velocity case.
We could try adding a velocity v again, but there’s an important catch to this. v must always be larger than vthr. Otherwise, it’s kinda strange to call this the high-speed regime when our boundary defining high speed is violated. Still, maybe we could do the following:
xp = 2 * xc * (ln(vi / vthr) - ln(v / vthr))
Amazingly, our guess again passes the x(500, 250) + x(250, 0) ==? x(500, 0) test. To be fair, I did try xp = 2 * xc * ln((vi - v) / vthr) as well, but this one didn’t pass the test at all. We’re on a roll again!
Our goals are of course to obtain t(x) and x(t) again. Only then can we estimate the depth of a bullet based on the remaining time step of the physics simulation tick, or the time until a certain distance (out of the obstacle?) is reached. Ummm. How do we do this? Maybe the time could again be easily described using t = avg(vi, v)/x? I tried, count with me!, the (1) arithmetic mean velocity, the (2) root-mean-squared velocity, the (3) geometric mean velocity, (4) inverse mean, the (5) logarithmic mean and basically all p-values of the (5 + infinity ~= infinity) stolarsky mean… And NONE were accurate (means listed in order, below).
avg(vi, v) =? (v,i + v) / 2
avg(vi, v) =? sqrt(v,i^2 / 2 + v^2 / 2)
avg(vi, v) =? sqrt(v,i * v)
avg(vi, v) =? 2 / (1 / v,i + 1 / v)
avg(vi, v) =? (v - v,i) / (ln(v) - ln(v,i))
avg(vi, v) =? ((v,i^p - v^p)/(p * v,i - v))^(1 / p - 1)
What I mean by none of them worked is that, when I try the tried- and tested t(500, 250) + t(250, 0) ==? t(500, 0) test, the results are never equal and each of these averaging methods fail my test!
Did I lose my mojo? Maybe, but I definitely lost my cool. Eventually, I decided to reach for a bottle of WolframAlpha… I mean, I entered the equation in WolframAlpha and observed the following:
xp = 2 * xc * (ln(vi / vthr) - ln(v / vthr))
ln(a) - ln(b) = ln(a/b)
xp = 2 * xc * ln((vi / vthr) / (v / vthr))
(a/b) / (c/b) = a/c
xp = 2 * xc * ln(vi / v)
vi / v = exp(xp / (2 * xc))
v = vi / exp(xp / (2 * xc)) = vi * exp(-xp / (2 * xc))
v = dx/dt = xp'
xp' = vi * exp(-xp / (2 * xc))
The WolframAlpha Gods had answered my prayers! Of course, they added a factor c1 again, of course. I immediately rewrote the entire piece of crap answer equation into t(xp).
xp(t) = 2 * xc * ln(-(c1 - t * vi) / (2 * xc))
xp / (2 * xc) = ln(-(c1 - t * vi) / (2 * xc))
-c1 + vi * t = 2 * xc * exp(xp / (2 * xc))
vi * t = 2 * xc * exp(xp / (2 * xc)) + c1
t = 2 * xc / vi * exp(xp / (2 * xc)) + c1 / vi
I couldn’t tell you why, but in a stroke of genius I remembered that units in an equation must match on both sides. So, t[s], meaning that c1 [?] / vi[m/s] = [? * s/m], it must be the case that c1 [m]. Secondly, I had already prepared all of these darned averages which weren’t accurate to more than up to 1%. Meaning, I could easily try a few variable combinations with meter-units until they more or less matched the other averaging methods! Wow. After about a minute, it was obvious that c1 [m] = -2 * xc.
xp(t) = 2 * xc * ln(-(-2 * xc - t * vi) / (2 * xc))
t(xp) = 2 * xc / vi * (exp(xp / (2 * xc)) - 1)
And that’s a wrap on this case! Well done everybody.
Non-deforming projectiles: Intermediate speed penetration
Penetration Depth x(t) [m] as Function of Time t [s] (all v)
x(t) = xc * ln((vi^2 + vthr^2) / (vthr^2) * (1 - 1 / (1 + tan(vthr * t / (2 * xc) + atan(vthr / vi)))^2))
Penetration Time t(x) [s] as Function of Depth x [m] (all v)
t(x) = 2 * xc / vthr * (atan(sqrt(1 / (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) - 1)) - atan(vthr / vi))
Maximum Distance xmax [m] and Maximum Time tmax [s] of Penetration
xmax = xc * ln(1 + vi^2 / (vthr^2))
tmax = xc / vthr * (pi() - 2 * atan(vthr/vi))
This derivation is the case of slight nightmares. The equations look ever so slightly familiar, but I can’t for the life of me derive this mess into a workable formula. Let’s do it anyways!
x = xc * ln(1 + vi^2 / vthr^2)
All variables are described before.. so let’s immediately get into adding a velocity v. Everything seems fine at this point! Gulp. All of the x(500, 250) + x(250, 0) ==? x(500, 0) tests work perfectly, so I guess we’re lucky with the guesses again.
xm = xc * (ln(1 + vi^2 / vthr^2) - ln(1 + v^2 / vthr^2))
Now to get the t(x) relation.. first off I hoped maybe one of the averages worked? And yea, no, no cigar. Similarly, the t(x) methods for low- and high-speed regimes were generally overshooting the value in t(500,250) + t(250,0) ==? t(500,0) tests. I looked to WolframAlpha for solutions, but the way I entered the formula changed the solutions I got – maybe I made an error in some of the operations deriving the differential equation. Therefore I ended up putting the whole thing in Mathematica, using DSolve[] and solving for x(t). Miraculously, one of the solutions was workable:
x(t) = xc * ln((vi^2 + vthr^2) * tan(vthr * (t + c1)/(2 * xc))^2 / (vthr^2 * (1 + tan(vthr * (t + c1)/(2 * xc))^2)
An absolutely terrible solution to the problem – too many characters! But, this form can be simplified a bit to make it more readable and workable:
x(t) = xc * ln((vi^2 + vthr^2) * tan(vthr * (t + c1)/(2 * xc))^2 / (vthr^2 * (1 + tan(vthr * (t + c1)/(2 * xc))^2)
x(t) = xc * ln((vi^2 + vthr^2) / (vthr^2) * f(t) / (1 + f(t))
f(t) / (1 + f(t) = ((1 + f(t)) - 1) / (1 + f(t) = 1 - 1 / (1 + f(t))
x(t) = xc * ln((vi^2 + vthr^2) / (vthr^2) * (1 - 1 / (1 + tan(vthr * (t + c1)/(2 * xc))^2))
Annoyingly, there’s another c1 [s] in this equation as well. We’ll leave it there for the time being, because it turns out to be quite a large formula. Instead, let us consider t(x), this is mainly a big rewriting job:
x(t) = xc * ln((vi^2 + vthr^2) / (vthr^2) * (1 - 1 / (1 + tan(vthr * (t + c1)/(2 * xc))^2)))
exp(x / xc) = (vi^2 + vthr^2) / (vthr^2) * (1 - 1 / (1 + tan(vthr * (t + c1)/(2 * xc))^2))
vthr^2 / (vi^2 + vthr^2) * exp(x / xc) = (1 - 1 / (1 + tan(vthr * (t + c1)/(2 * xc))^2))
- 1 / (1 + tan(vthr * (t + c1)/(2 * xc))^2) = vthr^2 / (vi^2 + vthr^2) * exp(x / xc) - 1
1 / (1 + tan(vthr * (t + c1)/(2 * xc))^2) = 1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)
1 + tan(vthr * (t + c1)/(2 * xc))^2 = 1 / (1 - vthr^2) / (vi^2 + vthr^2) * exp(x / xc))
tan(vthr * (t + c1)/(2 * xc))^2 = 1 / (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) - 1
tan(vthr * (t + c1)/(2 * xc)) = sqrt(1 / (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) - 1)
vthr * (t + c1)/(2 * xc) = atan(sqrt(1 / (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) - 1))
(t + c1) = 2 * xc / vthr * atan(sqrt(1 / (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) - 1))
t = 2 * xc / vthr * atan(sqrt(1 / (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) - 1)) - c1
t(x) = 2 * xc / vthr * atan(sqrt(1 / (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) - 1)) - c1
Now, the results from t(x) again seemed to be completely different from the other averages. Furthermore, if v == 0, the formula would return a division by zero-error. Looking into this specific error, it appeared that the contents of (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) were exactly zero when v = 0. Why? It seems that the following happened:
g(x(v)) = (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc))
x(v=0) = xc * ln(1 + vi^2 / vthr^2)
exp(xc * ln(1 + vi^2 / vthr^2) / xc) = exp(ln(1 + vi^2 / vthr^2)) = 1 + vi^2 / vthr^2
g(x(v=0)) = (1 - vthr^2 / (vi^2 + vthr^2) * (1 + vi^2 / vthr^2))
a / (b + a) * (1 + b / a) = a / (b + a) + a * b / (a * (b + a)) = a / (b + a) + b / (b + a)
= (a + b) / (b + a) = 1
g(x(v=0)) = (1 - 1) = 0
This issue can be circumvented. Whenever g(x(v~0)) ~= 0, 1/g(x) - 1 ~= Inf, such that sqrt(1/g(x) - 1) ~= Inf. Filling in atan(Inf) ~= Pi / 2. This then simplifies to t(x(v=0)) ~= xc / vthr * Pi - c1.
Another observation allowed me to find c1, namely that t(x) was again reversed, such that t(xmax) = 0 and t(0) = tmax. This again called for c1 = t(0), such that the following describes c1. The formula can be simplified a few times even:
c1 = t(0)
c1 = 2 * xc / vthr * atan(sqrt(1 / (1 - vthr^2 / (vi^2 + vthr^2)) - 1))
1 / (1 - a) - 1 = 1 / (1 - a) - (1 - a) / (1 - a) = (1 - 1 + a) / (1 - a) = a / (1 - a)
c1 = 2 * xc / vthr * atan(sqrt((vthr^2 / (vi^2 + vthr^2)) / (1 - vthr^2 / (vi^2 + vthr^2))))
(a / (a + b)) / (1 - a / (a + b)) = (a / (a + b)) / ((a + b) / (a + b) - a / (a + b)) = a / (a + b - a) = a / b
c1 = 2 * xc / vthr * atan(sqrt(vthr^2 / vi^2))
c1 = 2 * xc / vthr * atan(vthr / vi)
Filling in c1 at x(t) and t(x) then yields:
x(t) = xc * ln((vi^2 + vthr^2) / (vthr^2) * (1 - 1 / (1 + tan(vthr * (t + c1)/(2 * xc))^2))
tan(vthr * (t + c1)/(2 * xc))^2 = tan(vthr * (t + 2 * xc / vthr * atan(vthr / vi))/(2 * xc))^2
= tan(vthr * t / (2 * xc) + atan(vthr / vi)))^2
x(t) = xc * ln((vi^2 + vthr^2) / (vthr^2) * (1 - 1 / (1 + tan(vthr * t / (2 * xc) + atan(vthr / vi)))^2))
t(x) = 2 * xc / vthr * atan(sqrt(1 / (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) - 1)) - c1
t(x) = 2 * xc / vthr * atan(sqrt(1 / (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) - 1)) - 2 * xc / vthr * atan(vthr / vi)
t(x) = 2 * xc / vthr * (atan(sqrt(1 / (1 - vthr^2 / (vi^2 + vthr^2) * exp(x / xc)) - 1)) - atan(vthr / vi))
The maximum time and distance of the projectile can be described as:
xmax = xc * ln(1 + vi^2 / (vthr^2))
tmax = xc / vthr * (pi() - 2 * atan(vthr/vi))
Finally, the remaining velocity after distance x is described as:
v(x) = sqrt(vthr^2 * (exp(-x / xc) * (1 + vi^2 / vthr^2) - 1))
That was it
Go home. Have a bath. Get some rest. Here’s the fruits of our labour. As you can see, low speed (purple) and high speed (green) assumptions are no-good very-bad. On the other hand, the intermediate speed (pink) is gorgeous.

For the penetration timesteps, it really looks like the low speed results are fiiiiine while the high speed results quickly converge to the low speed ones.. and both converge to the intermediate speeds. Oh well, maybe it’s better to be certain.







