Simple cartridge model

Overview

The goal of this notebook is to create a simple model of firearm cartridge pressures and forces.

  • Simple Internal Ballistics Model
    • Chamber pressure curve
    • Bullet dynamics
    • Chamber wall friction
  • .22lr model
  • 9mm model
  • 357 model (TBD)
  • 40 S&W (TBD)

Simple Internal Ballistics

This simple model is a function of three main physical effects:

  1. Rapidly rising powder combustion pressure
  2. The combustion gas expanding as the bullet accelerates down the barrel
  3. The bolt thrust as a function of chamber pressure and case friction

Chamber Pressure Curve

The chamber pressure is governed by two main effects, the peak combustion pressure of the powder and the expansion of the gasses while the bullet accelerates down the barrel.

Chinn and Quickload tend to model the initial rise in powder combustion pressure as an exponential:

While Varmit Al shows some very nice pressure curves with more of an "S" shape to them initially:

For this model we choose to fit a Sigmoid Function:

$$S(t) = \frac{1}{1+\mathrm{e}^{-\alpha t}}$$

The curve below represents what the peak chamber pressure would be assuming that the bullet and bolt are frozen from moving. This is not the same as peak chamber pressure during normal operation. Once the bullet is allowed to accelerate down the barrel, the chamber pressures will be much smaller as the gas volume increases.


In [1]:
%pylab inline
from IPython.display import Math
from curves import sigmoid
import units as u

# sigmoid curve parameters
t_peak  = 5e-4                    # s, time for 50% sigmoid pressure
alpha_s = 6.0 / t_peak            # Combustion pressure timeconstant, 0-peak
p_peak  = 40 * u.ksi              # Peak combustion pressure

#Time axis
Ts = 1e-6
Tf = 3.0e-3
N = int(Tf / Ts)
t = array(range(0,N))*Ts

# Render the base curve
plot(t, p_peak/u.ksi*sigmoid(t, alpha_s, t_peak))
title('Combustion pressure vs time')
xlabel('Time, s')
ylabel('Pressure, ksi')
axis([0, 3e-3, 0, 45])
grid()
show()


Populating the interactive namespace from numpy and matplotlib

Bullet Dynamics

For the purposes of this analysis we will assume that the bullet is a mass pushed through a frictionless tube with a perfect gas seal. $$ F_{bullet} = P_{chamber} A_{bullet} $$ As the bullet accelerates down the barrel the pressure in the chamber will decrease as the gas expands. We will assume isentropic expansion which is governed by $$ \frac{p_1}{p_0} = \left( \frac{V_0}{V_1}\right) ^{\gamma} $$

Where $\gamma$ is the ratio of specific heats. Here $\gamma$ = 1.3 which should be in the ball park for hot air and $CO_2$.

22LR pressure curve

The following simulation combines the pressure curve with the bullet dynamics to create a more familiar chamber pressure cure. The model is calibrated to correlate with Varmit Al's 22lr data. The model assumes a 24.5" barrel, and a 40 grain .22lr bullet.


In [2]:
from cartridges import C22_Al
import units as u

c = C22_Al()

X = zeros([4,N])        # x, x_d, pressure, F_case
for i in range(N-1):
    x = X[0,i]
    x_d = X[1,i]
    x_dd = X[2,i] * c.a_chamber / c.m_bullet

    p_combust = c.p_peak*sigmoid(t[i], c.alpha_s, c.t_peak)
    if x >  24.5 * u.inch: 
        p_combust = 0
    p_k = p_combust*pow(c.v_case / (c.v_case + x*c.a_chamber), c.gama)

    F_case = p_k * c.a_chamber
    
    X[0,i+1] = x + x_d*Ts
    X[1,i+1] = x_d + x_dd*Ts
    X[2,i+1] = p_k
    X[3,i+1] = F_case

F_bolt = array(X[3,:])

figure(figsize=(12,8), dpi=150)
subplot(2,2,1)
plot(t,X[2,:] /u.ksi)
title('Chamber Pressure')
xlabel('Time, s')
ylabel('Chamber Pressure, ksi')
axis([0, 3e-3, 0, 8])
grid()

subplot(2,2,2)
plot(t, X[1,:]/u.ft)
title('Bullet Velocity')
xlabel('Time, s')
ylabel('Velocity, ft/s')
axis([0, 3e-3, 0, 1200])
grid()

show()


This model looks qualitatively pretty good. After choosing the rise time of the sigmoid function, the peak pressure was hand tuned until the the bullet velocity matched. Alternatively the pressures could have been tuned to match the peak pressure more closely at the expense of a lower bullet velocity.

Chamber wall pressure and friction

For bolt action and falling block actions, the above model is sufficient. For blowback or other actions where the case may move while under pressure, it is helpful to have an estimate for the friction between the case wall and the chamber wall. This model assumes that the chamber is smooth (not fluted.)

First we will look at the yield and burst pressures for the case using the thin wall hoop stress equation, then calculate the force between the brass and the chamber wall and finally estimate the frictional force acting between the case and the chamber.

The hoop stress equation is: $$ \sigma_{\theta} = \frac{Pd}{2t} $$ Where $\sigma_{\theta}$ is the stress on the brass, $P$ is the chamber pressure, $d$ is the chamber diameter and $t$ is the case wall thickness.

For brass, the yield stress is assumed to be ~29ksi and the ultimate strength is estimated to be ~80ksi. $$ P_{yield} = \frac{2\sigma_{yield}t}{d} \\ P_{burst} = \frac{2\sigma_{uts}t}{d} $$


In [3]:
p_yield = c.yield_case * 2.0*c.t_case / c.d_case
p_uts = c.uts_case * 2.0*c.t_case / c.d_case
print 'Case yield pressure: %.2f ksi' % (p_yield / u.ksi)
print 'Case burst pressure: %.2f ksi' % (p_uts / u.ksi)


Case yield pressure: 2.86 ksi
Case burst pressure: 7.86 ksi

Even for the relatively low velocity .22lr round, the case will deform significantly and rupture without support from the chamber.

The force acting between the case wall and the chamber wall is estimated as: $$ \begin{align} F_{wall} & = P_{case \_ wall} A_{case} \\ & = (P_{chamber} - P_{yield}) \pi d_{case} l_{case} \end{align} $$ Where $P_{case \_ wall}$ is the net pressure acting between the case and the chamber walls, $d_{case}$ is the diameter of the case and $l_{case}$ is the length of case acting on the wall (forming the cylindrical area).

The frictional force is then: $$ F_{friction} = \mu F_{wall} $$

We can now add this effect to the above simulation:


In [4]:
X = zeros([4,N])        # x, x_d, pressure, F_case_wall
for i in range(N-1):
    x = X[0,i]
    x_d = X[1,i]
    x_dd = X[2,i] * c.a_chamber / c.m_bullet

    p_combust = c.p_peak*sigmoid(t[i], c.alpha_s, c.t_peak)
    if x > 24.5 * u.inch:
        p_combust = 0
    p_k = p_combust*pow(c.v_case / (c.v_case + x*c.a_chamber), c.gama)

    F_case = p_k * c.a_chamber

    # Case wall friction removes force
    if p_k < p_yield:
        F_friction = 0
    else:
        a_case = pi*c.d_case*c.l_case
        F_friction = (p_k - p_yield) * a_case * c.u_cw

    F_case =p_k * c.a_chamber - F_friction
    if F_case < 0: 
        F_case = 0
    
    X[0,i+1] = x + x_d*Ts
    X[1,i+1] = x_d + x_dd*Ts
    X[2,i+1] = p_k
    X[3,i+1] = F_case

figure(figsize=(12,8), dpi=150)
subplot(2,2,1)
plot(t, F_bolt.T / u.lbf, label='w/o friction')
plot(t,X[3,:] / u.lbf, label='w friction')
legend(loc=0)
title('Bolt Force')
#xlabel('Time, s')
ylabel('Bolt Force, lbf')
axis([0, 2.5e-3, 0, 350])
grid()

subplot(2,2,2)
plot(t, X[1,:]/u.ft)
title('Bullet Velocity')
#xlabel('Time, s')
ylabel('Velocity, ft/s')
#axis([0, 2.5e-3, 0, 1200])
grid()

show()


The bolt force figure tells the story here. Once the pressure is high enough for the case to yield, the chamber begins supporting the brass and a frictional force develops between the two. Even with a relatively low coefficient of friction (0.1 in this case) the larger area of the case wall affords a relatively large force for the same chamber pressure when compared to the force on the case base. In this case of a .22lr the peak force applied to the bolt is reduced significantly by the friction between the case and the chamber.

9mm Cartridge Model

The process is repeated for the 9mm Luger cartridge below.


In [5]:
from cartridges import C9mm
import units as u

c = C9mm()

l_barrel = 16.0*u.inch
(t,Y) = c.sim(False, l_barrel)
(t,Z) = c.sim(True, l_barrel)

figure(figsize=(12,8), dpi=150)
subplot(2,2,1)
plot(Y[0,:] / u.inch ,Y[2,:] / u.ksi)
title('Chamber Pressure Vs Bullet Position')
ylabel('Pressure, kPsi')
axis([0, l_barrel / u.inch, 0, 40])
grid()

subplot(2,2,2)
plot(t, Y[2,:]/u.ksi)
title('Chamber Pressure vs Time')
ylabel('Pressure, kPsi')
axis([0, 1.5e-3, 0, 40])
grid()

subplot(2,2,3)
plot(Y[0,:] / u.inch, Y[1,:] / u.ft)
title('Bullet Velocity Vs Position')
xlabel('Bullet Position, inchs')
ylabel('Bullet Velocity, ft/s')
axis([0, l_barrel/u.inch, 0, 1600])
grid()

subplot(2,2,4)
plot(t, Y[1,:] / u.ft)
title('Bullet Velocity Vs Time')
xlabel('Time, s')
ylabel('Bullet Velocity ft/s')
axis([0, 1.5e-3, 0, 1600])
grid()

show()



In [6]:
figure(figsize=(12,8), dpi=150)
subplot(2,2,1)
plot(t,Y[3,:] / u.lbf, label='w/o friction')
plot(t,Z[3,:] / u.lbf, label='w friction')
legend(loc=0)
title('Bolt Force')
xlabel('Time, s')
ylabel('Bolt Force, lbf')
axis([0, 1.5e-3, 0, 4000])
grid()

subplot(2,2,2)
plot(t, Y[1,:]/u.ft)
title('Bullet Velocity')
xlabel('Time, s')
ylabel('Velocity, ft/s')
axis([0, 1.5e-3, 0, 1600])
grid()


Next steps

  • .38 Special and .357 magnum models
  • .40 S&W model

In [7]:
css_path = 'custom.css'
from ipynb_utils import inject_css2
from IPython.display import HTML
HTML(inject_css2(css_path))


Out[7]:

In [8]:
# %load_ext autoreload
# %autoreload 2
# %pylab inline