2.1 Interior ballistics

2.1.1 Control volume notation

m˙ijpsxiAigas
Figure 2.1: Control volume number i as represented in BlasterSim.

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 m˙ij, a projectile/plunger labeled p with mass mp,i, a spring with stiffness ki, mass ms,i, and spring precompression Δpre,i (see figure 1.4 for a graphical definition), cross-sectional area Ai, projectile/plunger position xi and velocity x˙i, and gas mass mi and energy Ei.

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, xi represents the projectile location, and ki and ms,i are set to zero;

  • a plunger tube for springers where p represents the plunger head, and xi represents the plunger head location;

  • a constant volume chamber for pneumatics where p represents a solid wall, mp,i is set to infinity, and ki and ms,i are set to zero;

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

  • a constant pressure atmosphere where xi is arbitrary, and ki and ms,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.4), whether a quantity is extensive or intensive can be determined from the units for that quantity, if unclear.

In the inputs, as diameters are provided, BlasterSim assumes that each control volume is cylindrical with a circular or annular cross-section. However, internally, BlasterSim uses only control volume cross-sectional area, not a diameter. BlasterSim internally makes no assumptions about the shape of the cross-section of each control volume, so BlasterSim can be extended to other control volume shapes in the future if desired, or an equivalent diameter with the intended cross-sectional area can be used in the meantime.

2.1.2 Conservation laws

BlasterSim internally tracks control volume mass and energy. The mass conservation equation for species k applied to a control volume i is [8, eq. 4.2, p. 147]:

dmk,idt=j(yk,jm˙jiyk,im˙ij), (2.1)

where mk,i is the mass of species k in control volume i, yk,i=mk,i/mi is the mass fraction of species k in control volume i, and m˙ji is the mass flow rate of all species from control volume j to control volume i.

The total mass in a control volume can be found by summing over the species:

mi=kmk,i (2.2)

The energy conservation equation for all gas species applied to a control volume i is [8, eq. 4.15, p. 157]:

dEg,idt=piAix˙i+j(m˙jihjm˙ijhi), (2.3)

where Eg,i is the total energy of the gas in control volume i, pi is the pressure in control volume i, Ai is the cross-sectional area of control volume i (barrel or plunger), x˙i is the velocity of the projectile/plunger in control volume i, and hi is the specific enthalpy of the gas in control volume i.

In the term m˙jihj, the use of specific enthalpy of the gas in another control volume (hj) comes from modeling the flow restriction as a “throttling device” which has negligible internal volume (removing the unsteady term), no internal work or heat transfer, and negligible gas kinetic and gravitational potential energy [8, p. 177]. For a throttling device, the specific enthalpy of the gas entering the flow restriction (hj) equals the specific enthalpy of the gas leaving the flow restriction, so the total enthalpy rate going from control volume j to control volume i is m˙jihj. The same logic applies for the m˙ijhi term.

Some simplifying assumptions are made in equation 2.3 that might be relaxed in the future. The gas is assumed to have no kinetic energy, so Eg,i contains only thermal energy and the kinetic energy flux term is set to zero. The gravitational potential energy of the gas is assumed to be zero, so the gravitational potential energy flux term is set to zero.

Given the lack of kinetic and gravitational potential energy, the total gas energy is related to only the internal energy in a control volume:

Eg,i=miui. (2.4)

The pressure work term [8, p. 39] in equation 2.3, piAix˙i intentionally has a leading negative sign. When a plunger is compressing the gas in the plunger tube, energy is entering the plunger tube gas. In this case, by BlasterSim convention, x˙i is negative, so the leading negative sign is canceled out, making the pressure work positive. Similarly, when the gas in the barrel is expanding, the pressure work should be negative, meaning energy is leaving the barrel gas. In this case, by BlasterSim convention, x˙i is positive, so the leading negative sign is kept, making the pressure work negative.

2.1.3 Equations of state

BlasterSim has two equations of state at present:

  • Ideal gas equation of state.

  • Constant specified pressure, temperature, and mass density equation of state.

The ideal gas equation of state gets the pressure from the following equation [8, eq. 3.33, p. 115]:

pi=miRiTiAixi. (2.5)

BlasterSim checks the validity of the ideal gas equation of state from the gas mixture critical pressure and will report an error if the ideal gas equation of state is used and invalid.

The constant specified pressure equation of state is used for the atmosphere, where the pressure, temperature, and mass density are expected to remain constant.

2.1.4 Thermodynamic properties

The specific heats cv,k and cp,k are assumed to be constant for each gas species. For each gas species, BlasterSim has the molar mass Mk, the ratio of specific heats γk, the specific internal energy u0,k, and specific enthalpy h0,k at reference temperature T0=300.0 K.

The species gas constant Rk is found from the molar mass Mk and universal gas constant R¯ via [8, p. 119]:

Rk=R¯Mk. (2.6)

From there, the specific heats at constant volume and constant pressure for species k can be found via [8, p. 119]:

cv,k =Rkγk1, (2.7)
cp,k =γkRkγk1. (2.8)

From there, the specific internal energy uk,i and specific enthalpy hk,i for a particular gas species k is found (via the constant specific heat assumption) through

uk,i(Ti) =u0,k+cv,k(TiT0,k), (2.9)
hk,i(Ti) =h0,k+cp,k(TiT0,k), (2.10)

where Ti is the temperature in control volume i.

For a particular control volume with a given array of mass fractions yk,i, the specific internal energy ui and specific enthalpy hi are found from the sums

ui(Ti) =kyk,iuk,i(Ti), (2.11)
hi(Ti) =kyk,ihk,i(Ti). (2.12)

2.1.5 Flow connection mass flow rate model

Control volumes can be connected so that mass can flow between them through flow restrictions. Flow restrictions include valves in the pneumatic case or any flow path with significant flow resistance. (If the flow resistance between two control volumes is negligible, it would make more sense to merge those control volumes into one.) BlasterSim uses a variation of the Beater [1, eq. 5.4, p. 46] compressible flow rate model modified to be smooth to support automatic differentiation (see § 2.2.2).

The mass flow rate from control volume i to control volume j is modeled as

m˙ij =αijm˙0,ij
+αijAe,ij(pigm˙(pjpi)pj)1bijRiTi
1fm˙(pjpi,bij), (2.13)

where αij is a term to account for valves opening (see § 2.1.6), m˙0 is a term used for testing only so that a specified mass flow rate can be used, Ae,ij is the effective area of the flow restriction, pi and pj are the pressures of control volumes i and j, respectively, gm˙ and fm˙ are functions defined below, bij is the critical pressure ratio of the flow restriction, Ri is the gas constant for the gas mixture in control volume i, and Ti is the temperature in control volume i.

For gm˙ and fm˙, let’s define pr=pj/pi and prl=0.999, where pr>prl is arbitrarily chosen to be where laminar flow occurs following the simplified laminar flow model of Beater [1, p. 46].

The purpose of gm˙ is to gradually and continuously turn on a term that is needed for pr near 1 to properly model laminar flow and also avoid unphysically high sensitivities to the pressure. The function gm˙ is

gm˙(pr)={0,prprl,2(1pr1prl)33(1pr1prl)2+1,pr>prl. (2.14)

The derivative of the pr>prl branch of gm˙ with respect to pr is zero at pr=prl so that the following function is continuous and continuously differentiable as it switched between the prprl and pr>prl branches.

fm˙ is intended to both account for laminar flow (pr>prl) and choked flow (pr<bij). First, a smooth minimum function must be defined [9]:

smin(x,y)=kln((exp(x/k)+exp(y/k))), (2.15)

where k=0.01, chosen arbitrarily. This is used as a continuous version of the conventional minimum function min(x,y).

Then fm˙ is

fm˙(pr,b)=(smin(pr,prl)b1b)212(1+tanh(prbps)), (2.16)

where ps=0.01, chosen arbitrarily to determine how quickly to switch between subsonic and choked flow.

2.1.6 Valve opening model

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. αij is the valve opening fraction, which is how open the valve is. αij=0 is fully closed m˙ij=0 and αij=1 is fully open. The specific equation BlasterSim uses is

αij(t)=α0,ij +α˙0,ij(ttopening,ij)
+(33α0,ij2α˙0,ij)(ttopening,ij)2
(22α0,ijα˙0,ij)(ttopening,ij)3, (2.17)

for t<topening,ij. For ttopening,ij, αij=1. α0,ij is the initial (time zero) valve opening fraction (useful when a valve isn’t used like in a springer). α˙0,ij is the initial valve opening rate. topening,ij is the valve opening time, the time it takes for the valve to fully open (reach αij=1).

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

αij(0) =α0,ij (set initial valve opening fraction) (2.18)
dαijdt(0) =α˙0,ij (set initial valve opening rate) (2.19)
αij(topening,ij) =1 (valve is fully open attopening,ij) (2.20)
dαijdt(topening,ij) =0 (needed for automatic differentiation) (2.21)

When the valve is fully open, αij no longer changes with time, so its derivative is zero. Consequently, the last constraint listed is needed to match the derivative at t=topening,ij, 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 (αij=α0,ij+(1α0,ij)t/topening,ij) from being used.

Note that to ensure monotonicity of αij (in other words, avoid oscillations of the valve opening fraction), α˙0,ij must satisfy the inequality 0α˙0,ij3(1α0,ij).

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

For pneumatics, α0=0 and α˙0=1. The second condition approximates a simple linear model. In the future α˙0 may become an input parameter if deemed necessary to improve accuracy. The SpudFiles Wiki [14] 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, α0=1 and α˙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.

2.1.7 Projectile and plunger equations of motion

The position xi of the projectile/plunger moves according to its velocity x˙i:

dxidt=x˙i. (2.22)

Newton’s second law applied to the projectile/plunger obtains its acceleration:

meff,idx˙idt=Ai(pipmirror,ipf,i)ki(xixmin,i+Δpre,i), (2.23)

where meff,i is the effective mass of the projectile/plunger in control volume i (defined below in this section, and set to be infinitely large for fixed volume control volumes so that the volume stays constant), Ai is the cross-sectional area of control volume i (barrel or plunger), pi is the pressure in control volume i, pmirror,i is the pressure in the control volume on the opposite side of the projectile/plunger of control volume i, pf,i is the pressure of friction (defined in § 2.1.8), ki is the spring stiffness of the plunger in control volume i (zero in control volumes where no spring is present like the barrel), xmin,i is the minimum location the projectile/plunger can take in control volume i (see § 2.1.9 for what happens when the projectile/plunger reaches this location), and Δpre,i is the precompression of the spring in control volume i (zero if no spring is present).

Gravitational force on the projectile/plunger is neglected for simplicity.

Note that to ease understanding of the CSV output, the variable xi in equation 2.23 in BlasterSim internally is different from x as printed by BlasterSim in the CSV output. See § 1.4 for details.

When a spring is present, the effective mass meff,i 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

meff,i=mp,i+Cmsms,i (2.24)

where mp,i is the mass of the plunger in control volume i, ms,i is the mass of the spring in control volume i, and Cms=13 as suggested by Ruby [11] for a stiff spring.

2.1.8 Projectile and plunger friction model

Consider a projectile/plunger p with both static and dynamic friction. BlasterSim uses an effective friction pressure pf,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 pfs,i and the magnitude of the dynamic friction pressure is pfd,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 dx˙idt=0; see equation 2.23), until the force exceeds the maximum static friction force. In other words, the actual static friction pressure pf0,i is the minimum of the equilibrium pressure pfe,i or the maximum static friction pressure, and in the direction of the equilibrium pressure. Mathematically, this can be stated as

pf0,i=sgn(pfe,i)min(pfs,i,|pfe,i|). (2.25)

The actual dynamic friction pressure is sgn(x˙i)pfd,i, taking into account both the magnitude of the dynamic friction pressure pfd,i and the direction.

The effective friction pressure pf,i then will switch between static and dynamic friction depending on the velocity:

pf,i={pf0,i,x˙i=0(static),sgn(x˙i)pfd,i,x˙i0(dynamic). (2.26)

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 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.26, it’s helpful to first rewrite equation 2.26 in terms of the sgn function:

pf,i=pf0,i+[pfd,isgn(x˙i)pf0,i]sgn(x˙i). (2.27)

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

sgn(x˙i)tanh(x˙ivs), (2.28)

where vs 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

pf,i=pf0,i+[pfd,itanh(x˙ivss)pf0,i]tanh(x˙ivsd). (2.29)

In equation 2.29, instead of a single vs, two different velocity scales vss and vsd are used to allow for finer control. Note that contrary to the non-smooth model (equation 2.26), the absolute value of equation 2.29 is not guaranteed to be at most pfs,i or pfd,i. To avoid the magnitude of the friction pressure from going above what a user likely intended, vss is made much smaller than vsd. 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 pf0,i is more difficult. Rather than replacing a single part of equation 2.25, an new function was constructed that does not resemble equation 2.25 algebraically but takes similar values. A piecewise function is implemented in cva.f90 as p_f0 so that the actual static friction pressure pf0,i equals the equilibrium pressure pfe,i for small values of pfe,i:

pf0,i={pfe,i,|pfe,i|ps,sgn(pfe,i)pf0,i,high,|pfe,i|>ps, (2.30)

where ps is chosen to be a somewhat arbitrary fraction of pfs,i and

pf0,i,high=pfs,itanh(|pfe,i|pspfs,ips+arctanh(pspfs,i)). (2.31)

pf0,i defined in this way will exactly balance static friction up to an equilibrium pressure of ps, and smoothly transition to the maximum static friction after that.

Note that despite the appearance of sgn(pfe,i) and |pfe,i|, equation 2.30 is continuous and has continuous derivatives as those terms only become discontinuous at pfe,i=0 and do not appear in the branch where pfe,i=0.

Equation 2.31 can be derived starting with a tanh function written in the following form:

pf0,i,high(pfe,i)=atanh(b|pfe,i|+c), (2.32)

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

pf0,i,high(ps) =ps (continuity of value) (2.33)
dpf0,i,highdpfe,i(ps) =1 (continuity of derivative) (2.34)
pf0,i,high() =pfs,i (maximum static friction pressure) (2.35)

2.1.9 Plunger impact

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, ep [2, pp. 524–525], which varies between 0 and 1. The velocity of the plunger before impact, x˙i, is related to the velocity of the plunger after impact, x˙i:

x˙i=ep,ix˙i. (2.36)

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, Δem,i is tracked in BlasterSim for each control volume:

Δem,i=12meff,i(x˙i2x˙i2). (2.37)

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

The location where plunger impact occurs is specified with xmin,i internally in BlasterSim, with xmin,i>0. By making xmin,i>0, BlasterSim avoids a problem where as xi goes to zero, mk,i and Ei also have to go to zero. That synchronization is difficult to enforce numerically without additional complications. The synchronization issue can be avoided by simply making xmin,i>0. In reality, all things have some dead volume, and where the plunger tube control volume ends and the barrel begins is arbitrary. In early versions of BlasterSim, all dead volume in springers was placed into the barrel, but xmin,i for the barrel and plunger tube are picked so that half the dead volume is in each.

Having a non-zero coefficient of restitution will lead to an infinite number of bounces with a decreasing time between each bounce. The number of bounces in reality is finite, and the decreasing time between bounces can cause similar numerical issues when the time between bounces becomes comparable to the time step. A simple and crude way to fix this problem is to set the velocity of the plunger after impact to zero if the magnitude of the velocity before impact is below a threshold. The threshold used in BlasterSim at the moment is 1.0 cm/s.

If the plunger velocity after impact is zero (either due to the coefficient of restitution being zero or the velocity threshold mentioned in the previous paragraph), BlasterSim can stall. After impact, the plunger 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. The only time step that would not push the plunger past the end of the plunger tube in this case is zero, causing simulation progress to stall because the same problem will happen at the next time iteration. BlasterSim at the moment does not apply a force to prevent the plunger from moving if it is at the end of the plunger tube, which would be more physically realistic. 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 real 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.