This simple model is a function of three main physical effects:
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()
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$.
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.
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)
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.
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()
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