2.1 Interior ballistics

Internally, BlasterSim simulate blasters using control volumes with arbitrary flow connections between control volumes. This allows BlasterSim to simulates spring and pneumatic blasters with the same core simulation code. This also allows simulating more atypical blasters without major code changes.

An isolated control volume numbered $i$ is shown in figure 2.1. A single flow connection to the control volume numbered $j$ is shown with mass flow rate ${\dot{m}}_{i\to j}$, a projectile/plunger labeled p with mass $m_{\text{p},i}$, a spring with stiffness $k_i$, mass $m_{\text{s},i}$, and spring precompression ${\mathrm{\Delta }}_{\text{pre},i}$ (see figure 1.2 for a graphical definition), cross-sectional area $A_i$, projectile/plunger position $x_i$ and velocity ${\dot{x}}_i$, and gas mass $m_i$ and energy $E_i$.

The integers $i$ and $j$ will be used as indices for control volumes.

Gas mass and energy as chosen as these are conserved variables and it is believed that using conserved variables will lead to better conservation properties in the numerical methods.

A control volume represented like this can be used to represent a variety of scenarios relevant to gas gun interior ballistics, including:

• •

  a barrel, where p represents the projectile, $x_i$ represents the projectile location, and $k_i$ and $m_{\text{s},i}$ are set to zero;

• •

  a plunger tube for springers where p represents the plunger head, and $x_i$ represents the plunger head location;

• •

  a constant volume chamber for pneumatics where p represents a solid wall, $m_{\text{p},i}$ is set to infinity, and $k_i$ and $m_{\text{s},i}$ are set to zero;

• •

  a piston tube used in HPA Nerf blasters similar to a plunger tube; and

• •

  a constant pressure atmosphere where $x_i$ is arbitrary, and $k_i$ and $m_{\text{s},i}$ are set to zero.

The gas mass could be composed of multiple species in BlasterSim. At the moment, the pneumatic and springer modes of BlasterSim use only dry air, but humidity and alternative gases are intended to be added in the future. For completeness, this guide will be written as if multiple gas species are present. The mass fraction $y$ and gas species index $k$ which identifies the gas are used.

In the typeset math in this documentation, when there are separate intensive and extensive versions of the same quantity, the two will be distinguished by making the extensive version uppercase and the intensive version lowercase. For example, internal energy per unit mass is $u$ and internal energy is $U$. Note that this convention is not followed in the typeset math when it is conventional to make an extensive quantity lowercase (like mass $m$) or an intensive quantity uppercase (like $R$), only when there are separate intensive and extensive versions of the same quantity. This capitalization convention is not followed in BlasterSim’s source code. See § 4.3 for a discussion of the capitalization convention followed in the source code. Because each variable in BlasterSim’s source code is given a unit (see § 3.1.3), whether a quantity is extensive or intensive can be determined from the units for that quantity, if unclear.

$$\frac{d{m}_{k,i}}{dt}=\sum _j\left(y_{k,j}{\dot{m}}_{j\to i}-y_{k,i}{\dot{m}}_{i\to j}\right)$$

(2.1)

$$\frac{d{E}_{\text{g},i}}{dt}=-p_i{A}_i{\dot{x}}_i+\sum _j\left({\dot{m}}_{j\to i}{h}_j-{\dot{m}}_{i\to j}{h}_i\right)$$

(2.2)

$$m_i=\sum _k{m}_{k,i}$$

(2.3)

$$E_i=m_i{u}_i$$

(2.4)

$$p_i=\frac{m_i{R}_i{T}_i}{A_i{x}_i}$$

(2.5)

$$u_{k,i}(T)=u_{0,k}+c_{\text{v},k}(T_i-T_{0,k})$$

(2.6)

$$h_{k,i}(T)=h_{0,k}+c_{\text{p},k}(T_i-T_{0,k})$$

(2.7)

Beater [1, ch. 5]

Valves do not open instantaneously, and sometimes slow opening has a significant impact on performance. BlasterSim models how far a valve is open with a prescribed valve opening model. ${\alpha }_{i\to j}$ is the valve opening fraction, which is how open the valve is. ${\alpha }_{i\to j}=0$ is fully closed ${\dot{m}}_{i\to j}=0$ and ${\alpha }_{i\to j}=1$ is fully open. The specific equation BlasterSim uses is

	${\alpha }_{i\to j}={\alpha }_{0,i\to j}$	$+{\dot{\alpha }}_{0,i\to j}\left({\displaystyle \frac{t}{t_{\text{opening},i\to j}}}\right)$
		$+(3-3{\alpha }_{0,i\to j}-2{\dot{\alpha }}_{0,i\to j}){\left({\displaystyle \frac{t}{t_{\text{opening},i\to j}}}\right)}^2$
		$-(2-2{\alpha }_{0,i\to j}-{\dot{\alpha }}_{0,i\to j}){\left({\displaystyle \frac{t}{t_{\text{opening},i\to j}}}\right)}^3,$		(2.8)

for . For $t\ge t_{\text{opening},i\to j}$, ${\alpha }_{i\to j}=1$. ${\alpha }_{0,i\to j}$ is the initial (time zero) valve opening fraction (useful when a valve isn’t used like in a springer). ${\dot{\alpha }}_{0,i\to j}$ is the initial valve opening rate. $t_{\text{opening},i\to j}$ is the valve opening time, the time it takes for the valve to fully open (reach ${\alpha }_{i\to j}=1$).

This equation may seem overly complicated, but is the simplest polynomial equation that satisfies a few constraints:

	${\alpha }_{i\to j}(0)$	$={\alpha }_{0,i\to j}$	(set initial valve opening fraction)		(2.9)
	${\displaystyle \frac{d{\alpha }_{i\to j}}{dt}}(0)$	$={\dot{\alpha }}_{0,i\to j}$	(set initial valve opening rate)		(2.10)
	${\alpha }_{i\to j}(t_{\text{opening},i\to j})$	$=1$	$\text{(valve is fully open at}{t}_{\text{opening},i\to j}\text{)}$		(2.11)
	${\displaystyle \frac{d{\alpha }_{i\to j}}{dt}}(t_{\text{opening},i\to j})$	$=0$	(needed for automatic differentiation)		(2.12)

When the valve is fully open, ${\alpha }_{i\to j}$ no longer changes with time, so its derivative is zero. Consequently, the last constraint listed is needed to match the derivative at $t=t_{\text{opening},i\to j}$, which is necessary for automatic differentiation. See § 2.2.2 for more about automatic differentiation in BlasterSim. The last constraint prevents a simple linear model (${\alpha }_{i\to j}={\alpha }_{0,i\to j}+(1-{\alpha }_{0,i\to j})t/t_{\text{opening},i\to j}$) from being used.

Note that to ensure monotonicity of ${\alpha }_{i\to j}$ (in other words, avoid oscillations of the valve opening fraction), ${\dot{\alpha }}_{0,i\to j}$ must satisfy the inequality $0\le {\dot{\alpha }}_{0,i\to j}\le 3(1-{\alpha }_{0,i\to j})$.

The pneumatic and springer cases will now be discussed, with the $i\to j$ subscript dropped for simplicity as there is only one flow restriction in both cases.

For pneumatics, ${\alpha }_0=0$ and ${\dot{\alpha }}_0=1$. The second condition approximates a simple linear model. In the future ${\dot{\alpha }}_0$ may become an input parameter if deemed necessary to improve accuracy. The SpudFiles Wiki [11] gives some estimated opening times:

• •

  Burst disks: likely under 1 ms

• •

  Pilot-operated valves (like QEVs, “back-pressure tanks”, “cores”): 3–5 ms

• •

  Ball valves: about 100 ms if hand activated

These values are recommended as starting points only. For modeling any particular blaster, it is better to try to independently determine the valve opening time through something like high speed video.

For springers, ${\alpha }_0=1$ and ${\dot{\alpha }}_0=0$. As there is no valve in a springer, this simply sets the flow restriction to always be open.

This simple valve opening model is most accurate for manually operated valves. More detailed modeling of pilot-operated valves could make determining the valve opening time unnecessary, but this alternative valve opening model has not yet been added to BlasterSim.

	${\displaystyle \frac{d{x}_i}{dt}}$	$={\dot{x}}_i$		(2.13)
	${\displaystyle \frac{d{\dot{x}}_i}{dt}}$	$={\displaystyle \frac{A_i}{m_{\text{eff},i}}}\left(p_i-p_{\text{mirror},i}-p_{\text{f},i}\right)-{\displaystyle \frac{k_i}{m_{\text{eff},i}}}\left(x_i-x_{\text{min},i}+{\mathrm{\Delta }}_{\text{pre},i}\right)$		(2.14)

When a spring is present, the effective mass of the plunger factors in the spring mass. The spring is not moving at a uniform velocity as one end is stationary, so it would be incorrect to add all the spring mass to the effective mass. The effective mass equation used is

$$m_{\text{eff},i}=m_{\text{p},i}+C_{\text{ms}}{m}_{\text{s},i}$$

(2.15)

where $C_{\text{ms}}=\frac{1}{3}$ as suggested by Ruby [9] for a stiff spring.

Note that the variable $x_i$ in equation 2.14 internally in BlasterSim differs from x as printed by BlasterSim in the CSV output. To ease understanding for users of BlasterSim, x in the CSV output is defined as the distance between the projectile/plunger location and the minimum $x$ location like the end of the plunger tube. Mathematically, $\text{x}=x_i-x_{\text{min},i}$.

Consider a projectile/plunger p with both static and dynamic friction. BlasterSim uses an effective friction pressure $p_{\text{f},i}$ defined to oppose the direction of motion. Using a pressure instead of force allows for more familiar units and is conventional in spud guns. The maximum static friction pressure is $p_{\text{fs},i}$ and the magnitude of the dynamic friction pressure is $p_{\text{fd},i}$. When static friction is active, the actual static friction force is just enough force to equilibrate the forces so that p doesn’t move (where $\frac{d{\dot{x}}_i}{dt}=0$; see equation 2.14), until the force exceeds the maximum static friction force. In other words, the actual static friction pressure $p_{\text{f0},i}$ is the minimum of the equilibrium pressure $p_{\text{fe},i}$ or the maximum static friction pressure, and in the direction of the equilibrium pressure. Mathematically, this can be stated as

$$p_{\text{f0},i}=\mathrm{sgn}(p_{\text{fe},i})\cdot \min (p_{\text{fs},i},|p_{\text{fe},i}|).$$

(2.16)

The actual dynamic friction pressure is $\mathrm{sgn}({\dot{x}}_i)\cdot p_{\text{fd},i}$, taking into account both the magnitude of the dynamic friction pressure $p_{\text{fd},i}$ and the direction.

The effective friction pressure $p_{\text{f},i}$ then will switch between static and dynamic friction depending on the velocity:

(2.17)

BlasterSim does not use this model because it is not smooth. Smooth functions are useful for automatic differentiation (see § 2.2.2) and beneficial for numerical stability. Also, the use of strict floating point equality in the ${\dot{x}}_i=0$ branch can lead to numerical reproducibility issues. BlasterSim uses a smoothed version of this model which avoids these problems.

To create a smoothed version of equation 2.17, it’s helpful to first rewrite equation 2.17 in terms of the $\mathrm{sgn}$ function:

$$p_{\text{f},i}=p_{\text{f0},i}+[p_{\text{fd},i}-\mathrm{sgn}({\dot{x}}_i)\cdot p_{\text{f0},i}]\cdot \mathrm{sgn}({\dot{x}}_i).$$

(2.18)

The smooth approximation to equation 2.17 will be made by noting that

$$\mathrm{sgn}({\dot{x}}_i)\approx \tanh \left(\frac{{\dot{x}}_i}{v_{\text{s}}}\right),$$

(2.19)

where $v_{\text{s}}$ is a small velocity scale over which the $\tanh$ function changes from -1 to 1.

Consequently, the function implemented in cva.f90 as p_f is

$$p_{\text{f},i}=p_{\text{f0},i}+\left[p_{\text{fd},i}-\tanh \left(\frac{{\dot{x}}_i}{v_{\text{ss}}}\right)\cdot p_{\text{f0},i}\right]\cdot \tanh \left(\frac{{\dot{x}}_i}{v_{\text{sd}}}\right).$$

(2.20)

In equation 2.20, instead of a single $v_{\text{s}}$, two different velocity scales $v_{\text{ss}}$ and $v_{\text{sd}}$ are used to allow for finer control. Note that contrary to the non-smooth model (equation 2.17), the absolute value of equation 2.20 is not guaranteed to be at most $p_{\text{fs},i}$ or $p_{\text{fd},i}$. To avoid the magnitude of the friction pressure from going above what a user likely intended, $v_{\text{ss}}$ is made much smaller than $v_{\text{sd}}$. The bounds violation comes from both $\tanh$ terms changing at the same time, so it is believed making the two scales very different so that one finishes changing sooner than the other will avoid this problem.

Making a smooth approximation for $p_{\text{f0},i}$ is more difficult. Rather than replacing a single part of equation 2.16, an new function was constructed that does not resemble equation 2.16 algebraically but takes similar values. A piecewise function is implemented in cva.f90 as p_f0 so that the actual static friction pressure $p_{\text{f0},i}$ equals the equilibrium pressure $p_{\text{fe},i}$ for small values of $p_{\text{fe},i}$:

(2.21)

where $p_{\text{s}}$ is chosen to be a somewhat arbitrary fraction of $p_{\text{fs},i}$ and

$$p_{\text{f0},i,\text{high}}=p_{\text{fs},i}\tanh \left(\frac{|p_{\text{fe},i}|-p_{\text{s}}}{p_{\text{fs},i}-p_{\text{s}}}+\mathrm{arctanh}\left(\frac{p_{\text{s}}}{p_{\text{fs},i}}\right)\right).$$

(2.22)

$p_{\text{f0},i}$ defined in this way will exactly balance static friction up to an equilibrium pressure of $p_{\text{s}}$, and smoothly transition to the maximum static friction after that.

Note that despite the appearance of $\mathrm{sgn}(p_{\text{fe},i})$ and $|p_{\text{fe},i}|$, equation 2.21 is continuous and has continuous derivatives as those terms only become discontinuous at $p_{\text{fe},i}=0$ and do not appear in the branch where $p_{\text{fe},i}=0$.

Equation 2.22 can be derived starting with a $\tanh$ function written in the following form:

$$p_{\text{f0},i,\text{high}}(p_{\text{fe},i})=a\tanh \left(b\cdot |p_{\text{fe},i}|+c\right),$$

(2.23)

and finding the values of the coefficients $a$, $b$, and $c$ that satisfy the following constraints:

	$p_{\text{f0},i,\text{high}}(p_{\text{s}})$	$=p_{\text{s}}$	(continuity of value)		(2.24)
	${\displaystyle \frac{d{p}_{\text{f0},i,\text{high}}}{d{p}_{\text{fe},i}}}(p_{\text{s}})$	$=1$	(continuity of derivative)		(2.25)
	$p_{\text{f0},i,\text{high}}(\mathrm{\infty })$	$=p_{\text{fs},i}$	(maximum static friction pressure)		(2.26)

When the plunger impacts the end of the plunger tube, the plunger tube will rebound. The amount of rebound can be controlled through the coefficient of restitution, $e_{\text{p}}$ [2, pp. 524–525], which varies between 0 and 1. The velocity of the plunger before impact, ${\dot{x}}_i$, is related to the velocity of the plunger after impact, ${\dot{x}}_i^{\prime }$:

$${\dot{x}}_i^{\prime }=-e_{\text{p},i}{\dot{x}}_i.$$

(2.27)

The negative sign shows that the plunger will change from moving towards the end of the plunger tube to moving away from the end of the plunger tube.

The amount of energy dissipated in the impact, $\mathrm{\Delta }{e}_{\text{m},i}$ is tracked in BlasterSim for each control volume:

$$\mathrm{\Delta }{e}_{\text{m},i}=\frac{1}{2}{m}_{\text{eff},i}\left({\dot{x}}_i^2-{\dot{x}}_i^{\prime 2}\right).$$

(2.28)

Tracking the total amount of energy dissipated in impact can be useful to improve blaster longevity.

If the coefficient of restitution is zero so that the plunger velocity after impact is zero, BlasterSim applies an approximation for the moment. The plunger will become immobile after impact. This approximation is a simple way to prevent simulation progress from stalling. After impact, the plunger tube is located at the end of the plunger tube. If the coefficient of restitution is zero, then the plunger will initially have zero velocity after impact. However, the force on the plunger will push it past the end of the plunger tube in the simulation. BlasterSim simply identifies when the plunger moves past the end of the plunger tube and picks a time step to exactly hit the end of the plunger tube. BlasterSim at the moment does not additionally apply an extra force to prevent the plunger from moving if it is at the end of the plunger tube, which would be more physically realistic. Consequently, to prevent the plunger tube from moving, it is made immobile in this situation. In the future, if desirable, BlasterSim could apply a force when the plunger is at the end of the plunger tube to keep the plunger stationary. The difference with the force is that it might be possible for the force to be overcome by the pressure inside the plunger, making the plunger start moving again. Making the plunger immobile prevents that motion, but is a reasonable approximation for springers. It would not be a reasonable approximation when simulating the full cycle of a HPA blaster using an air cylinder. It is likely that coefficients of restitution that are too small may have similar stalling problems that this approximation will not handle properly.