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 is shown in figure 2.1. A single flow connection to the control volume numbered is shown with mass flow rate , a projectile/plunger labeled p with mass , a spring with stiffness , mass , and spring precompression (see figure 1.4 for a graphical definition), cross-sectional area , projectile/plunger position and velocity , and gas mass and energy .
The integers and 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, represents the projectile location, and and are set to zero;
a plunger tube for springers where p represents the plunger head, and represents the plunger head location;
a constant volume chamber for pneumatics where p represents a solid wall, is set to infinity, and and are set to zero;
a piston tube used in HPA Nerf blasters similar to a plunger tube; and
a constant pressure atmosphere where is arbitrary, and and 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 and gas species index 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 and internal energy is . Note that this convention is not followed in the typeset math when it is conventional to make an extensive quantity lowercase (like mass ) or an intensive quantity uppercase (like ), 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.
BlasterSim internally tracks control volume mass and energy. The mass conservation equation for species applied to a control volume is [8, eq. 4.2, p. 147]:
| (2.1) |
where is the mass of species in control volume , is the mass fraction of species in control volume , and is the mass flow rate of all species from control volume to control volume .
The total mass in a control volume can be found by summing over the species:
| (2.2) |
The energy conservation equation for all gas species applied to a control volume is [8, eq. 4.15, p. 157]:
| (2.3) |
where is the total energy of the gas in control volume , is the pressure in control volume , is the cross-sectional area of control volume (barrel or plunger), is the velocity of the projectile/plunger in control volume , and is the specific enthalpy of the gas in control volume .
In the term , the use of specific enthalpy of the gas in another control volume () 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 () equals the specific enthalpy of the gas leaving the flow restriction, so the total enthalpy rate going from control volume to control volume is . The same logic applies for the 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 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:
| (2.4) |
The pressure work term [8, p. 39] in equation 2.3, 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, 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, is positive, so the leading negative sign is kept, making the pressure work negative.
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]:
| (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.
The specific heats and are assumed to be constant for each gas species. For each gas species, BlasterSim has the molar mass , the ratio of specific heats , the specific internal energy , and specific enthalpy at reference temperature K.
The species gas constant is found from the molar mass and universal gas constant via [8, p. 119]:
| (2.6) |
From there, the specific heats at constant volume and constant pressure for species can be found via [8, p. 119]:
| (2.7) | ||||
| (2.8) |
From there, the specific internal energy and specific enthalpy for a particular gas species is found (via the constant specific heat assumption) through
| (2.9) | ||||
| (2.10) |
where is the temperature in control volume .
For a particular control volume with a given array of mass fractions , the specific internal energy and specific enthalpy are found from the sums
| (2.11) | ||||
| (2.12) |
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 to control volume is modeled as
| (2.13) |
where is a term to account for valves opening (see § 2.1.6), is a term used for testing only so that a specified mass flow rate can be used, is the effective area of the flow restriction, and are the pressures of control volumes and , respectively, and are functions defined below, is the critical pressure ratio of the flow restriction, is the gas constant for the gas mixture in control volume , and is the temperature in control volume .
For and , let’s define and , where is arbitrarily chosen to be where laminar flow occurs following the simplified laminar flow model of Beater [1, p. 46].
The purpose of is to gradually and continuously turn on a term that is needed for near 1 to properly model laminar flow and also avoid unphysically high sensitivities to the pressure. The function is
| (2.14) |
The derivative of the branch of with respect to is zero at so that the following function is continuous and continuously differentiable as it switched between the and branches.
is intended to both account for laminar flow () and choked flow (). First, a smooth minimum function must be defined [9]:
| (2.15) |
where , chosen arbitrarily. This is used as a continuous version of the conventional minimum function .
Then is
| (2.16) |
where , chosen arbitrarily to determine how quickly to switch between subsonic and choked flow.
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. is the valve opening fraction, which is how open the valve is. is fully closed and is fully open. The specific equation BlasterSim uses is
| (2.17) |
for . For , . is the initial (time zero) valve opening fraction (useful when a valve isn’t used like in a springer). is the initial valve opening rate. is the valve opening time, the time it takes for the valve to fully open (reach ).
This equation may seem overly complicated, but is the simplest polynomial equation that satisfies a few constraints:
| (set initial valve opening fraction) | (2.18) | ||||
| (set initial valve opening rate) | (2.19) | ||||
| (2.20) | |||||
| (needed for automatic differentiation) | (2.21) |
When the valve is fully open, no longer changes with time, so its derivative is zero. Consequently, the last constraint listed is needed to match the derivative at , 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 () from being used.
Note that to ensure monotonicity of (in other words, avoid oscillations of the valve opening fraction), must satisfy the inequality .
The pneumatic and springer cases will now be discussed, with the subscript dropped for simplicity as there is only one flow restriction in both cases.
For pneumatics, and . The second condition approximates a simple linear model. In the future 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, and . 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.
The position of the projectile/plunger moves according to its velocity :
| (2.22) |
Newton’s second law applied to the projectile/plunger obtains its acceleration:
| (2.23) |
where is the effective mass of the projectile/plunger in control volume (defined below in this section, and set to be infinitely large for fixed volume control volumes so that the volume stays constant), is the cross-sectional area of control volume (barrel or plunger), is the pressure in control volume , is the pressure in the control volume on the opposite side of the projectile/plunger of control volume , is the pressure of friction (defined in § 2.1.8), is the spring stiffness of the plunger in control volume (zero in control volumes where no spring is present like the barrel), is the minimum location the projectile/plunger can take in control volume (see § 2.1.9 for what happens when the projectile/plunger reaches this location), and is the precompression of the spring in control volume (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 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 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
| (2.24) |
where is the mass of the plunger in control volume , is the mass of the spring in control volume , and as suggested by Ruby [11] for a stiff spring.
Consider a projectile/plunger p with both static and dynamic friction. BlasterSim uses an effective friction pressure 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 and the magnitude of the dynamic friction pressure is . 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 ; see equation 2.23), until the force exceeds the maximum static friction force. In other words, the actual static friction pressure is the minimum of the equilibrium pressure or the maximum static friction pressure, and in the direction of the equilibrium pressure. Mathematically, this can be stated as
| (2.25) |
The actual dynamic friction pressure is , taking into account both the magnitude of the dynamic friction pressure and the direction.
The effective friction pressure then will switch between static and dynamic friction depending on the velocity:
| (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 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 function:
| (2.27) |
The smooth approximation to equation 2.26 will be made by noting that
| (2.28) |
where is a small velocity scale over which the function changes from -1 to 1.
Consequently, the function implemented in cva.f90 as p_f is
| (2.29) |
In equation 2.29, instead of a single , two different velocity scales and 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 or . To avoid the magnitude of the friction pressure from going above what a user likely intended, is made much smaller than . The bounds violation comes from both 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 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 equals the equilibrium pressure for small values of :
| (2.30) |
where is chosen to be a somewhat arbitrary fraction of and
| (2.31) |
defined in this way will exactly balance static friction up to an equilibrium pressure of , and smoothly transition to the maximum static friction after that.
Note that despite the appearance of and , equation 2.30 is continuous and has continuous derivatives as those terms only become discontinuous at and do not appear in the branch where .
Equation 2.31 can be derived starting with a function written in the following form:
| (2.32) |
and finding the values of the coefficients , , and that satisfy the following constraints:
| (continuity of value) | (2.33) | ||||
| (continuity of derivative) | (2.34) | ||||
| (maximum static friction pressure) | (2.35) |
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, [2, pp. 524–525], which varies between 0 and 1. The velocity of the plunger before impact, , is related to the velocity of the plunger after impact, :
| (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, is tracked in BlasterSim for each control volume:
| (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 internally in BlasterSim, with . By making , BlasterSim avoids a problem where as goes to zero, and 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 . 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 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 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.