1 Introduction

Over the years, the approach to solving the main problem of internal ballistics has evolved significantly [1,2,3,4,5,6,7]. Nowadays, most design processes are computer-aided, but in the case of internal ballistics they are not generally common available. The analytical solution of the system of nonlinear differential equations describing the phases of the phenomenon of a shot in a gun barrel does not exist. Therefore, in order to obtain more accurate solutions to the system of equations, it is worth using numerical methods [8,9,10,11,12,13,14,15]. On the basis of the adopted physical model, a mathematical model is created, which should take into account the specific features of a given problem. The solution to the main problem of the internal ballistics of barrel weapons for the adopted physical model consists in determining the following parameters: total pressure, dynamic pressure, static pressure, average static pressure, temperature of gaseous products of gunpowder combustion (GCP), velocity of the projectile inside the barrel, displacement of the projectile inside the barrel.

The solution of the internal ballistics problem has been divided into main four phases, each of which is characterized by its own assumptions and is different from the others: preliminary phase, initial phase, pyrodynamic phase, adiabatic phase. Detailed descriptions of these phases are presented in detail by authors of the publications [3, 6, 15,16,17].

The ballistic phases presented above are described by differential equations and analytical equations. In the next section of the publication, the system of equations for individual phases is presented, which together build a mathematical model for solving the main problem of internal ballistics in the numerical approach proposed in the further part of this publication.

2 Mathematical model [18, 19]

The list of basic physical quantities used in the model:

p–pressure,

v–velocity,

l–displacement,

\(\rho _{m}\)–bulk density of GPM,

\(l_{0} \)–length of the cartridge chamber,

\(W_{0}\)–volume of the cartridge chamber.

\(W_{s}\)–volume of free space behind the bottom of the projectile,

\(L\left( t \right) =l_{0}+l\left( t \right) =\frac{W_{0}}{s}+l\left( t \right) \),

\(u_{1}\left[ \frac{m}{sPa} \right] \)–proper burn rate,

\(\varphi \)–projectile mass fictitiousness coefficient,

s–cross-sectional area of the barrel.

\(p_\textrm{s}\left( x, t \right) \)–static pressure depend of x from the range \(0\le x\le L\left( t \right) \) at the given t,

\(p_\textrm{d}\left( x, t \right) \)–dynamic pressure,

\(p_\textrm{a}\)–ambient pressure,

\(p_{z}\)–ignition pressure of the propellant charge,

\(t_{z}\)–time of ignition of the propellant charge,

\(S\left( e \right) \left[ m^{2} \right] \)–function of the combustion surface of powder grains,

\(e_{s}\left( t \right) \left[ m \right] \)–the thickness of the burnt powder grain layer depends on time.

\(S_{1}\left[ m^{2} \right] \) oraz \(\mathrm {\Lambda }_{\textrm{1}} \left[ m^{3} \right] \)–the surface area and the initial volume of the powder grain,

\(\omega \)–weight of the powder charge,

\(\omega _{z}\)–weight of igniter,

\(t_{z_{k}}\)–time of igniters decomposition into gaseous combustion products,

\(\delta \)–bulk powder mass density.

\(\varphi \)–fictitious mass of the projectile,

\(G_{Z}=\frac{\omega _{z}}{t_{z_{k}}}\)–the average mass flow of gaseous combustion products of the igniter,

G–the average mass flow of gaseous combustion products of the gunpowder,

\(\psi \)–the relative mass of the gunpowder burnt

\(\psi _{Z}\)–the relative mass of the igniter burnt,

\({\dot{N}}\)–the rate of creation of relative masses of gaseous combustion products,

\(\mathrm {\Lambda }_{\textrm{s}}\)–initial powder charge volume,

\(\alpha \)–covolume of the GCP,

R–gas constant,

T–temperature.

\(q_{\vartheta }\)–explosion heat of the propellant,

\(q_{z}\)–explosion heat of the igniter.

2.1 Preliminary phase

Preliminary phase that lasts in the pressure region:

$$\begin{aligned} p_{a}\le p\left( t \right) \le p_{z} \end{aligned}$$

and the corresponding time interval:

$$\begin{aligned} 0\le t\le t_\textrm{zap} \end{aligned}$$

is described by the following set of equations:

$$\begin{aligned}{} & {} \frac{\textrm{d}e_{s}\left( t \right) }{\textrm{d}t}={\dot{e}}_{s}\left( t \right) =0 \end{aligned}$$
(2.1)
$$\begin{aligned}{} & {} \frac{\textrm{d}\mathrm {\Lambda }_{\textrm{s}}\left( t \right) }{\textrm{d}t}=\dot{\Lambda }_{s}\left( t \right) =0 \end{aligned}$$
(2.2)
$$\begin{aligned}{} & {} \frac{\textrm{d}\psi \left( t \right) }{\textrm{d}t}={\dot{\psi }}\left( t \right) =0 \end{aligned}$$
(2.3)
$$\begin{aligned}{} & {} \frac{\textrm{d}\psi _{z}\left( t \right) }{\textrm{d}t}={\dot{\psi }}_{z}\left( t \right) =\frac{G_{z}}{\omega _{z}} \end{aligned}$$
(2.4)
$$\begin{aligned}{} & {} \frac{\textrm{d}N\left( t \right) }{\textrm{d}t}={\dot{N}}\left( t \right) =\dot{\psi _{z}}{\left( t \right) \omega }_{z} \end{aligned}$$
(2.5)
$$\begin{aligned}{} & {} \frac{\textrm{d}v\left( t \right) }{\textrm{d}t}={\dot{v}}\left( t \right) =0 \end{aligned}$$
(2.6)
$$\begin{aligned}{} & {} \frac{\textrm{d}l\left( t \right) }{\textrm{d}t}={\dot{l}}\left( t \right) =0 \end{aligned}$$
(2.7)
$$\begin{aligned}{} & {} \frac{\textrm{d}T}{\textrm{d}t}={\dot{T}}=\frac{k-1}{R}\left( \frac{p_{s}\left( L, t \right) }{{{(\rho }_{m}\left( t \right) )}^{2}}\frac{\textrm{d}\rho _{m}\left( t \right) }{\textrm{d}t}+q_{\vartheta }\frac{\textrm{d}\psi \left( t \right) }{\textrm{d}t}+q_{z}\frac{d\psi _{z}\left( t \right) }{\textrm{d}t} \right) \end{aligned}$$
(2.8)
$$\begin{aligned}{} & {} \rho _{m}\left( t \right) =0 \end{aligned}$$
(2.9)
$$\begin{aligned}{} & {} W_{s}\left( t \right) =W_{o}-\alpha \omega _{z}\psi _{z}\left( t \right) \end{aligned}$$
(2.10)
$$\begin{aligned}{} & {} p\left( t \right) =R\frac{N\left( t \right) T\left( t \right) }{W_{s}\left( t \right) } \end{aligned}$$
(2.11)
$$\begin{aligned}{} & {} \rho \left( t \right) =\frac{\omega _{s}^{*}}{W_{s}\left( t \right) } \end{aligned}$$
(2.12)
$$\begin{aligned}{} & {} p_{d}=\frac{\rho _{m}\left( v\left( t \right) \right) ^{2}}{2} \end{aligned}$$
(2.13)
$$\begin{aligned}{} & {} p_{s}\left( L, t \right) =p\left( t \right) -\frac{\rho _{m}\left( v\left( t \right) \right) ^{2}}{2} \end{aligned}$$
(2.14)
$$\begin{aligned}{} & {} {<p_{s}\left( x, t \right) >}_{L\left( t \right) }=p\left( t \right) -\frac{\rho _{m}\left( t \right) }{6}\left( v\left( t \right) \right) ^{2} \end{aligned}$$
(2.15)

where: \(\mathrm {\Lambda }_{\textrm{s}}\)—initial powder charge volume, R—gas constant of GCP, \(\omega ^{*}=\omega +\omega _{z}\), \(L\left( t \right) =l_{o}+l\)

2.2 Initial phase

Initial phase that lasts in the pressure region:

$$\begin{aligned} p_{z}\le p\le p_{0} \end{aligned}$$

and the corresponding time interval:

$$\begin{aligned} t_{z}\le t\le t_{0} \end{aligned}$$

is described by the following set of equations:

$$\begin{aligned}{} & {} \frac{\textrm{d}e_{s}\left( t \right) }{\textrm{d}t}={\dot{e}}_{s}\left( t \right) =u_{1}p\left( t \right) \end{aligned}$$
(2.16)
$$\begin{aligned}{} & {} \frac{\textrm{d}\mathrm {\Lambda }_{\textrm{s}}\left( t \right) }{\textrm{d}t}=\dot{\Lambda }_{s}(t)=S(e)\frac{de_{s}\left( t \right) }{\textrm{d}t}[\frac{m^{3}}{s}] \end{aligned}$$
(2.17)
$$\begin{aligned}{} & {} \frac{\textrm{d}\psi \left( t \right) }{\textrm{d}t}={\dot{\psi }}\left( t \right) =\frac{\delta }{\omega }S\left( e \right) \frac{de_{s}\left( t \right) }{dt}=\frac{\delta }{\omega }{\dot{\Lambda }}_{s}\left( t \right) \end{aligned}$$
(2.18)
$$\begin{aligned}{} & {} \frac{\textrm{d}\psi _{z}\left( t \right) }{\textrm{d}t}={\dot{\psi }}_{z}\left( t \right) =\frac{G_{z}}{\omega _{z}} \left[ \frac{1}{s} \right] \end{aligned}$$
(2.19)
$$\begin{aligned}{} & {} \frac{\textrm{d}N\left( t \right) }{\textrm{d}t}={\dot{N}} \left( t \right) =\dot{\psi _{z}}{\left( t \right) \omega }_{z}+\dot{\psi }\left( t \right) \omega \end{aligned}$$
(2.20)
$$\begin{aligned}{} & {} \frac{\textrm{d}v\left( t \right) }{\textrm{d}t}={\dot{v}}\left( t \right) =0 \end{aligned}$$
(2.21)
$$\begin{aligned}{} & {} \frac{\textrm{d}l\left( t \right) }{\textrm{d}t}={\dot{l}}\left( t \right) =0 \end{aligned}$$
(2.22)
$$\begin{aligned}{} & {} \frac{\textrm{d}T}{\textrm{d}t}={\dot{T}}=\frac{k-1}{R}\left( \frac{p_{s}\left( L, t \right) }{{{(\rho }_{m}\left( t \right) )}^{2}}\frac{\textrm{d}\rho _{m}\left( t \right) }{\textrm{d}t}+q_{\vartheta }\frac{\textrm{d}\psi \left( t \right) }{\textrm{d}t}+q_{z}\frac{d\psi _{z}\left( t \right) }{\textrm{d}t} \right) \end{aligned}$$
(2.23)
$$\begin{aligned}{} & {} \rho _{m}\left( t \right) =0 \end{aligned}$$
(2.24)
$$\begin{aligned}{} & {} W_{s}\left( t \right) =W_{o}-\frac{\omega }{\mathrm {\delta }}\left( 1-\psi \left( t \right) \right) -\alpha \omega \psi \left( t \right) -\alpha \omega _{z}\psi _{z}\left( t \right) \end{aligned}$$
(2.25)
$$\begin{aligned}{} & {} p\left( t \right) =R\frac{N\left( t \right) T\left( t \right) }{W_{s}\left( t \right) } \end{aligned}$$
(2.26)
$$\begin{aligned}{} & {} \rho \left( t \right) =\frac{\omega _{s}^{*}}{W_{s}\left( t \right) } \end{aligned}$$
(2.27)
$$\begin{aligned}{} & {} p_{d}=\frac{\rho _{m}\left( v\left( t \right) \right) ^{2}}{2} \end{aligned}$$
(2.28)
$$\begin{aligned}{} & {} p_{s}\left( L, t \right) =p\left( t \right) -\frac{\rho _{m}\left( v\left( t \right) \right) ^{2}}{2} \end{aligned}$$
(2.29)
$$\begin{aligned}{} & {} {<p_{s}\left( x, t \right) >}_{L\left( t \right) }=p\left( t \right) -\frac{\rho _{m}\left( t \right) }{6}\left( v\left( t \right) \right) ^{2} \end{aligned}$$
(2.30)

and after reaching \(t_{z_{k}}\):

$$\begin{aligned}{} & {} \frac{\textrm{d}e_{s}\left( t \right) }{\textrm{d}t}={\dot{e}}_{s}\left( t \right) =u_{1}p\left( t \right) \end{aligned}$$
(2.31)
$$\begin{aligned}{} & {} \frac{\textrm{d}\mathrm {\Lambda }_{\textrm{s}}\left( t \right) }{\textrm{d}t}=\dot{\Lambda }_{s}(t)=S(e)\frac{\textrm{d}e_{s}\left( t \right) }{dt} \end{aligned}$$
(2.32)
$$\begin{aligned}{} & {} \frac{\textrm{d}\psi \left( t \right) }{\textrm{d}t}={\dot{\psi }}\left( t \right) =\frac{\delta }{\omega }S\left( e \right) \frac{\textrm{d}e_{s}\left( t \right) }{\textrm{d}t}=\frac{\delta }{\omega }{\dot{\Lambda }}_{s}\left( t \right) \end{aligned}$$
(2.33)
$$\begin{aligned}{} & {} \frac{\textrm{d}\psi _{z}\left( t \right) }{\textrm{d}t}={\dot{\psi }}_{z}\left( t \right) =0 \end{aligned}$$
(2.34)
$$\begin{aligned}{} & {} \frac{\textrm{d}N\left( t \right) }{\textrm{d}t}={\dot{N}}\left( t \right) =\dot{\psi _{z}}{\left( t \right) \omega }_{z}+{\dot{\psi }}\left( t \right) \omega \end{aligned}$$
(2.35)
$$\begin{aligned}{} & {} \frac{\textrm{d}v\left( t \right) }{\textrm{d}t}={\dot{v}}\left( t \right) =0 \end{aligned}$$
(2.36)
$$\begin{aligned}{} & {} \frac{\textrm{d}l\left( t \right) }{\textrm{d}t}={\dot{l}}\left( t \right) =0 \end{aligned}$$
(2.37)
$$\begin{aligned}{} & {} \frac{\textrm{d}T}{\textrm{d}t}={\dot{T}}=\frac{k-1}{R}\left( \frac{p_{s}\left( L, t \right) }{{{(\rho }_{m}\left( t \right) )}^{2}}\frac{\textrm{d}\rho _{m}\left( t \right) }{\textrm{d}t}+q_{\vartheta }\frac{\textrm{d}\psi \left( t \right) }{\textrm{d}t}+q_{z}\frac{\textrm{d}\psi _{z}\left( t \right) }{\textrm{d}t} \right) \end{aligned}$$
(2.38)
$$\begin{aligned}{} & {} \rho _{m}\left( t \right) =0 \end{aligned}$$
(2.39)
$$\begin{aligned}{} & {} W_{s}\left( t \right) =W_{o}+sl\left( t \right) -\frac{\omega }{\mathrm {\delta }}\left( 1-\psi \left( t \right) \right) -\alpha \omega \psi \left( t \right) -\alpha \omega _{z}\psi _{z}\left( t \right) \end{aligned}$$
(2.40)
$$\begin{aligned}{} & {} p\left( t \right) =R\frac{N\left( t \right) T\left( t \right) }{W_{s}\left( t \right) } \end{aligned}$$
(2.41)
$$\begin{aligned}{} & {} \rho \left( t \right) =\frac{\omega _{s}^{*}}{W_{s}\left( t \right) } \end{aligned}$$
(2.42)
$$\begin{aligned}{} & {} p_{\textrm{d}}=\frac{\rho _{m}\left( v\left( t \right) \right) ^{2}}{2} \end{aligned}$$
(2.43)
$$\begin{aligned}{} & {} p_{s}\left( L, t \right) =p\left( t \right) -\frac{\rho _{m}\left( v\left( t \right) \right) ^{2}}{2} \end{aligned}$$
(2.44)
$$\begin{aligned}{} & {} {<p_{s}\left( x, t \right) >}_{L\left( t \right) }=p\left( t \right) -\frac{\rho _{m}\left( t \right) }{6}\left( v\left( t \right) \right) ^{2} \end{aligned}$$
(2.45)

2.3 Pyrodynamic phase

Pyrodynamic phase is described by the range of the layer thickness of the burnt powder grain at time \(t_{o}\) (i.e., the time when the gaseous combustion products-GCP reach the forcing pressure), and the thickness of the burnt grain layer at time \(t_{k}\) (the time when the grain burn finish):

$$\begin{aligned} e_{s}\left( t_{o} \right) \le e_{s}\left( t \right) \le e_{s}\left( t_{k} \right) \end{aligned}$$

The system of equations for this phase will take the form:

$$\begin{aligned}{} & {} \frac{\textrm{d}e_{s}\left( t \right) }{\textrm{d}t}={\dot{e}}_{s}(t)=u_{1}{<p_{s}\left( x, t \right) >}_{L\left( t \right) } \end{aligned}$$
(2.46)
$$\begin{aligned}{} & {} \frac{\textrm{d}\mathrm {\Lambda }_{\textrm{s}}\left( t \right) }{\textrm{d}t}=\dot{\Lambda }_{s}(t)=S(e)\frac{\textrm{d}e_{s}\left( t \right) }{\textrm{d}t}[\frac{m^{3}}{s}] \end{aligned}$$
(2.47)
$$\begin{aligned}{} & {} \frac{\textrm{d}\psi \left( t \right) }{\textrm{d}t}={\dot{\psi }}\left( t \right) =\frac{\delta }{\omega }S\left( e \right) \frac{\textrm{d}e_{s}\left( t \right) }{\textrm{d}t}=\frac{\delta }{\omega }{\dot{\Lambda }}_{s}\left( t \right) \end{aligned}$$
(2.48)
$$\begin{aligned}{} & {} \frac{\textrm{d}\psi _{z}\left( t \right) }{\textrm{d}t}={\dot{\psi }}_{z}\left( t \right) =0 \end{aligned}$$
(2.49)
$$\begin{aligned}{} & {} \frac{\textrm{d}N\left( t \right) }{\textrm{d}t}={\dot{N}}\left( t \right) =\dot{\psi _{z}}{\left( t \right) \omega }_{z}+{\dot{\psi }}\left( t \right) \omega \end{aligned}$$
(2.50)
$$\begin{aligned}{} & {} \frac{\textrm{d}v\left( t \right) }{\textrm{d}t}={\dot{v}}\left( t \right) =\frac{Sp_{s}\left( L, t \right) }{\varphi m} \end{aligned}$$
(2.51)
$$\begin{aligned}{} & {} \frac{\textrm{d}l\left( t \right) }{\textrm{d}t}={\dot{l}}\left( t \right) =v\left( t \right) \end{aligned}$$
(2.52)
$$\begin{aligned}{} & {} \frac{\textrm{d}T}{\textrm{d}t}={\dot{T}}=\frac{k-1}{R}\left( \frac{p_{s}\left( L, t \right) }{{{(\rho }_{m}\left( t \right) )}^{2}}\frac{\textrm{d}\rho _{m}\left( t \right) }{\textrm{d}t}+q_{\vartheta }\frac{\textrm{d}\psi \left( t \right) }{\textrm{d}t}+q_{z}\frac{\textrm{d}\psi _{z}\left( t \right) }{\textrm{d}t} \right) \end{aligned}$$
(2.53)
$$\begin{aligned}{} & {} \rho _{m}\left( t \right) =\frac{\omega ^{*}}{W\left( t \right) }=\frac{\omega ^{*}}{sL\left( t \right) } \end{aligned}$$
(2.54)
$$\begin{aligned}{} & {} W_{s}\left( t \right) =W_{o}+sl\left( t \right) -\frac{\omega }{\mathrm {\delta }}\left( 1-\psi \left( t \right) \right) -\alpha \omega \psi \left( t \right) -\alpha \omega _{z}\psi _{z}\left( t \right) \end{aligned}$$
(2.55)
$$\begin{aligned}{} & {} p\left( t \right) =R\frac{N\left( t \right) T\left( t \right) }{W_{s}\left( t \right) } \end{aligned}$$
(2.56)
$$\begin{aligned}{} & {} \rho \left( t \right) =\frac{\omega _{s}^{*}}{W_{s}\left( t \right) } \end{aligned}$$
(2.57)
$$\begin{aligned}{} & {} p_{d}=\frac{\rho _{m}\left( v\left( t \right) \right) ^{2}}{2} \end{aligned}$$
(2.58)
$$\begin{aligned}{} & {} p_{s}\left( L, t \right) =p\left( t \right) -\frac{\rho _{m}\left( v\left( t \right) \right) ^{2}}{2} \end{aligned}$$
(2.59)
$$\begin{aligned}{} & {} {<p_{s}\left( x, t \right) >}_{L\left( t \right) }=p\left( t \right) -\frac{\rho _{m}\left( t \right) }{6}\left( v\left( t \right) \right) ^{2} \end{aligned}$$
(2.60)

2.4 Adiabatic phase

Adiabatic phase is limited by the variability of the projectile displacement in the gun tube:

$$\begin{aligned} l\left( t_{k} \right) \le l\left( t \right) \le l_{w} \end{aligned}$$

where: \(l\left( t_{k} \right) \)—location of the projectiles bottom at the time of the end of burning the powder grains \(t_{k}\), \(l_{w}\)—gun tube length.

The system of equations for this phase will take the form:

$$\begin{aligned}{} & {} \frac{\textrm{d}e_{s}\left( t \right) }{\textrm{d}t}={\dot{e}}_{s}\left( t \right) = 0 \end{aligned}$$
(2.61)
$$\begin{aligned}{} & {} \frac{\textrm{d}\mathrm {\Lambda }_{\textrm{s}}\left( t \right) }{\textrm{d}t}={\dot{\Lambda }}_{s}\left( t \right) = 0 \end{aligned}$$
(2.62)
$$\begin{aligned}{} & {} \frac{\textrm{d}\psi \left( t \right) }{\textrm{d}t}={\dot{\psi }}\left( t \right) =0 \end{aligned}$$
(2.63)
$$\begin{aligned}{} & {} \frac{\textrm{d}\psi _{z}\left( t \right) }{\textrm{d}t}={\dot{\psi }}_{z}\left( t \right) =0 \end{aligned}$$
(2.64)
$$\begin{aligned}{} & {} \frac{\textrm{d}N\left( t \right) }{\textrm{d}t}={\dot{N}}\left( t \right) =0 \end{aligned}$$
(2.65)
$$\begin{aligned}{} & {} \frac{\textrm{d}v\left( t \right) }{\textrm{d}t}={\dot{v}}\left( t \right) =\frac{Sp_{s}\left( L, t \right) }{\varphi m} \end{aligned}$$
(2.66)
$$\begin{aligned}{} & {} \frac{\textrm{d}l\left( t \right) }{\textrm{d}t}={\dot{l}}\left( t \right) =v\left( t \right) \end{aligned}$$
(2.67)
$$\begin{aligned}{} & {} \frac{\textrm{d}T}{\textrm{d}t}={\dot{T}}=\frac{k-1}{R}\left( \frac{p_{s}\left( L, t \right) }{{{(\rho }_{m}\left( t \right) )}^{2}}\frac{\textrm{d}\rho _{m}\left( t \right) }{\textrm{d}t}+q_{\vartheta }\frac{\textrm{d}\psi \left( t \right) }{\textrm{d}t}+q_{z}\frac{\textrm{d}\psi _{z}\left( t \right) }{\textrm{d}t} \right) \end{aligned}$$
(2.68)
$$\begin{aligned}{} & {} \rho _{m}\left( t \right) =\frac{\omega ^{*}}{W\left( t \right) }=\frac{\omega ^{*}}{sL\left( t \right) } \end{aligned}$$
(2.69)
$$\begin{aligned}{} & {} W_{s}\left( t \right) =W_{o}+sl\left( t \right) -\frac{\omega }{\mathrm {\delta }}\left( 1-\psi \left( t \right) \right) -\alpha \omega \psi \left( t \right) -\alpha \omega _{z}\psi _{z}\left( t \right) \end{aligned}$$
(2.70)
$$\begin{aligned}{} & {} p\left( t \right) =R\frac{N\left( t \right) T\left( t \right) }{W_{s}\left( t \right) } \end{aligned}$$
(2.71)
$$\begin{aligned}{} & {} \rho \left( t \right) =\frac{\omega _{s}^{*}}{W_{s}\left( t \right) } \end{aligned}$$
(2.72)
$$\begin{aligned}{} & {} p_{d}=\frac{\rho _{m}\left( v\left( t \right) \right) ^{2}}{2} \end{aligned}$$
(2.73)
$$\begin{aligned}{} & {} p_{s}\left( L, t \right) =p\left( t \right) -\frac{\rho _{m}\left( v\left( t \right) \right) ^{2}}{2} \end{aligned}$$
(2.74)
$$\begin{aligned}{} & {} {<p_{s}\left( x, t \right) >}_{L\left( t \right) }=p\left( t \right) -\frac{\rho _{m}\left( t \right) }{6}\left( v\left( t \right) \right) ^{2} \end{aligned}$$
(2.75)

3 Program for numerical solution of internal ballistics

Using the equations formulated in the previous section describing the mathematical solution to the main problem of internal ballistics in Lagrange coordinates and using the Delphi 7 programming environment [20], the “Ball-Lagr” numerical program was created. It enables numerical simulation of the internal ballistics of classic barrel weapons. Delphi is described as a universal programming language and software supplies an integrated development environment. The Delphi dialect based on the Object Pascal programming language. In “Ball-Lagr,” analytic equations have been defined as functions that are executed in each iteration of program execution. The differential equations, on the other hand, are divided into functions in individual ballistic phases and only during a given phase is the part of the equations corresponding to it performed. Operation algorithm of the Ball-Lagr program is shown on the following schema (Fig. 1).

Fig. 1
figure 1

Block diagram of the general algorithm of the Ball-Lagr program

Full size image

Digital tool allows users to generate the graphs of each of the 15 equations, especially the pressure versus time course and the projectile velocity versus time. From the programming point of view, eight differential equations are kept in the state vector, when the remaining equations using procedures to perform calculations while solutions to differential equations are generating.

Numerical simulations determine changes in values over time of the following ballistic parameters such as:

  • Burnt grain layer density \(e_{s}\),

  • Volume of burned grains \(\mathrm {\Lambda }_{\textrm{s}}\),

  • The relative mass of burnt powder grains \(\psi \) and igniter \(\psi _{z}\),

  • Mass GCP \(\omega _{s}\),

  • Velocity v and displacement l of the projectiles bottom in the gun tube,

  • Volumetric density \(\rho _{m}\) GPM oraz \(\rho \) GCP,

  • Free volume \(W_{s}\) behind the bottom of the projectile and the volume of burnt grains \(\mathrm {\Lambda }_{\textrm{s}}\),

  • Temperature T of the GCP,

  • Total pressure at the bottom of the cartridge chamber p,

  • Dynamic pressure \(p_{d}\) MGP at the bottom of the projectile and the displacement l of the projectiles bottom in the gun tube,

  • Static pressure \(p_{s}\) GPM on the projectiles bottom and projectile velocity v,

  • Average static pressure \(p_{s}\) of GCP behind projectiles bottom.

Launching the Ball-Lagr program requires specifying the parameters controlling the program and the initial values for 8 differential equations, which are specified in the “Oblicz” (Calculate) table, as well as data on the physico-chemical and ballistic parameters of the propellants (BMP), shown in Figs. 2 and 3.

Fig. 2
figure 2

Ball-Lagr user window for entering control parameters and initial values

Full size image
Fig. 3
figure 3

Ball-Lagr user window for entering powder data

Full size image

where:

\(e_{1}\)—the thickness of the combustible layer of powder,

\(l_{w}\)— the path traveled by the bottom of the projectile in the gun tube,

eps, eps1, eps2—a variable for the precision function, whose role was to compact the step in a certain area,

skok—step for the RK4 (Runge–Kutta methods) calculation procedure between successive iterations of the program,

\(p_{z}\)—ignition pressure value,

\(t_\textrm{zap}\)—time to reach ignition pressure,

a1—variable for selecting the combustion surface (grain shape). Takes values from the range \(< 0\div 3>\). For the variable value \(a_{1}:= 0\), the "PowSpal" function calculates the burning area of spherical powder grains. For \(a_{1}:=1\), cylinder-shaped combustion surface of grains. Dla \(a_{1}:=2\) tube-shaped combustion surface. For \(a_{1}:=3\) or powder grains in the shape of cuboids,

wek.x—vector of initial conditions for differential equations.

where:

\(G_{z}\)—average mass flow of gaseous combustion products of the igniter

\(\alpha \)—covolume,

k—isentrope exponent,

\(q_{z}\)—heat associated with the igniter,

R—gas constant,

\(p_{0}\)—forcing pressure,

\(T_{o}\)—ambient temperature,

\(u_{1}\)—proper burn rate,

\(\varphi \)—secondary work factor, factor of the fictitious mass of the projectile,

\(\mathrm {\delta }\)—bulk powder mass density,

\(m_{p}\)—mass of projectile,

d—caliber,

\(p_{a}\)—ambient pressure,

\(q_{p}\)— heat associated with the powder,

\(\omega \)—mass of the powder.

Parameters related to grain dimensions (data downloaded depending on the adopted value of the a\(_{\textrm{1}}\) variable from the calculations tab):

\(D_{z}\)—grain diameter (in the case of the spherical-shape grain variant),

\(d_{w}\)—cylindrical diameter (in the case of tube-shaped grain variant),

\(l_{ziar}\)—grain length,

a, b, c—cuboid dimensions (in the case of cuboid-shaped grain variant).

\(\omega _{z}\)—mass of igniter,

\(l_{o}\)—chamber length.

Having the ballistic parameters of the selected type of powder, as well as the tool which is the designed Ball-Lagr program, it was possible to carry out a simulation for the selected powder weight and the selected type of weapon, which allowed to obtain internal ballistics solution diagrams. The tool provides a wide range of possibilities for analyzing the entered data. In the next chapter compares the results of shooting tests and numerical simulations for the 7.62x39 mm small arms cartridge.

4 Pyrodynamic research and simulation evaluation

The purpose of this study was to experimentally obtain the course of pressure as a function of time and the muzzle velocity of the bullet when firing cartridges of 7.62x39mm caliber from different manufacturers, with variable powder charge weights, and to compare them with the results of numerical simulations performed using the Ball-Lagr program.

The measuring equipment consisted of: the UPB universal ballistic stand with the AK-47 rifle, a set of control and measurement apparatus for measuring the pressure of powder gases using the piezo-electric method (Kistler), a set of control and measurement apparatus for measuring the projectile velocity (photoelectric gates).

Figure (Fig. 4) shows the arrangement of the pressure sensors that recorded the pressure waveforms at the given measuring points, i.e., sensor no. 1—pressure at the cartridge chamber, sensor no. 2—in the barrel duct, sensor no. 3—at the gas chamber. As part of this work, the measurement from sensor No. 1 was used for the analysis.

Fig. 4
figure 4

The research stand with an AK-47 rifle and indicated pressure measurement sensors

Full size image

The experimentally obtained measurements will be presented in the form of maximum pressure values for individual tests together with the muzzle velocity of the projectile. The results will be presented for individual types of powder, taking into account the weight of the powder sample. Next, a digital simulation of the shot phenomenon will be presented, which will be analyzed with experimental results.

The criterion for evaluating the quality of the digital solution to the main problem of the internal ballistics of classic barreled weapons in the paper was a comparison of the statistical results of experimental studies of the solution to the main problem of the internal ballistics of classic barreled weapons at selected points characteristic of the course of pressure and projectile velocity, with characteristic points on the curves of pressure and projectile velocity development for digital solution.

4.1 Results of ballistic measurements for SM ammunition cal. \(7.62 \times 39\) mm FMJ standard charge

Table 1 presents the results of ballistic measurements for a standard mass of the powder charge for a series of 5 shots

Table 1 The results of ballistic measurements for a standard mass of the powder charge
Full size table

The following figures show an experimental diagram of the pressure of the powder gases (GSP) in the cartridge chamber (Fig. 5) and the analogous graph of the pressure of the powder gases and the velocity of the projectile in the barrel obtained from numerical simulations (Fig. 6) for chosen shot.

Fig. 5
figure 5

An experimental diagram of the pressure of the powder gases (GSP) in the cartridge chamber (standard charge)

Full size image
Fig. 6
figure 6

Numerical simulation of the pressure and bullet velocity in the barrel duct as a function of time (standard charge)

Full size image

Based on the obtained simulation results and the results of ballistic measurements, it can be seen that the maximum value of the total pressure at the bottom of the chamber p obtained by simulation is 292.413 MPa and is within a range of three standard deviations from the mean value, i.e., 289.612 ± 9.066 MPa. The time to reach the maximum pressure in the simulation was 1.423 ms, while the time to reach maximum pressure based on ballistic tests is in the range 1.440 ± 0.017 ms. The value of the muzzle velocity obtained from the simulation is 663.085 m/s and is within a band of three standard deviations from the average value of the experiment, i.e., 666.018 ± 4.199 m/s. The shape of the graph is similar to the one obtained from the experiment. Summing up, the simulation gives a satisfactory reflection of the results of the internal ballistics solution by the method of experimental research.

4.2 Results of ballistic measurements for SM ammunition cal. \(7.62 \times 39\) mm FMJ reduced charge

Table 2 presents the results of ballistic measurements for a reduced mass of the powder charge for a series of 5 shots.

Table 2 The results of ballistic measurements for a reduced mass of the powder charge
Full size table

The following figures show an experimental diagram of the pressure of the powder gases (GSP) in the cartridge chamber (Fig. 7) and the analogous graph of the pressure of the powder gases and the velocity of the projectile in the barrel obtained from numerical simulations (Fig. 8) for chosen shot (reduced charge).

Fig. 7
figure 7

An experimental diagram of the pressure of the powder gases (GSP) in the cartridge chamber (reduced charge)

Full size image
Fig. 8
figure 8

Numerical simulation of the pressure and bullet velocity in the barrel duct as a function of time (reduced charge)

Full size image

Based on the obtained simulation results and the results of ballistic measurements, it can be seen that the maximum value of the total pressure at the bottom of the chamber p obtained by simulation is 151.650 MPa and is within a range of three standard deviations from the mean value, i.e., 133.398 ± 23.114 MPa. The time to reach the maximum pressure in the simulation was 1.510 ms, while the time to reach maximum pressure based on ballistic tests is in the range 1.537 ± 0.044 ms. The value of the muzzle velocity obtained from the simulation is 482.608 m/s and is within a band of three standard deviations from the average value of the experiment, i.e., 489.490 ± 22.497 m/s. Analogous to the previous results, the shape of the graph is similar to the one obtained from the experiment. Summing up, the simulation gives a satisfactory reflection of the results of the internal ballistics solution by the method of experimental research.

4.3 Results of ballistic measurements for SM ammunition cal. \(7.62 \times 39\) mm FMJ increased charge

Table 3 presents the results of ballistic measurements for a increased mass of the powder charge for a series of 5 shots.

Table 3 The results of ballistic measurements for a increased mass of the powder charge
Full size table
Fig. 9
figure 9

An experimental diagram of the pressure of the powder gases (GSP) in the cartridge chamber (increased charge)

Full size image
Fig. 10
figure 10

Numerical simulation of the pressure and bullet velocity in the barrel duct as a function of time (increased charge)

Full size image

The following figures show an experimental diagram of the pressure of the powder gases (GSP) in the cartridge chamber (Fig. 9) and the analogous graph of the pressure of the powder gases and the velocity of the projectile in the barrel obtained from numerical simulations (Fig. 10) for chosen shot (increased charge).

Based on the obtained simulation results and the results of ballistic measurements, it can be seen that the maximum value of the total pressure at the bottom of the chamber p obtained by simulation is 362.416 MPa and is within a range of three standard deviations from the mean value, i.e., 330.664 ± 40.611 MPa. The time to reach the maximum pressure in the simulation was 1.421 ms, while the time to reach maximum pressure based on ballistic tests is in the range 1.428 ± 0.007 ms. The value of the muzzle velocity obtained from the simulation is 700.201 m/s and is within a band of three standard deviations from the average value of the experiment, i.e., 693.260 ± 16.840 m/s. Analogous to the previous results, the shape of the graph is similar to the one obtained from the experiment. Summing up, the simulation gives a satisfactory reflection of the results of the internal ballistics solution by the method of experimental research.

5 Conclusions

Conclusions that can be drawn on the basis of performed numerical and experimental analyses are as follows:

  1. 1.

    Based on the presented mathematical model, a digital tool was created to solve the main problem of internal ballistics, which was designed in the Delphi 7 programming environment.

  2. 2.

    The digital tool, based on the entered input data (gun powder and projectile data), generates a waveform of a given quantity determined on the basis of each of the 15 equations as a function of time and allows you to export this data to a spreadsheet.

  3. 3.

    As a result of the experimental tests carried out for the \(7.62 \times 39\) mm FMJ cartridge, i.e., the measurement of the pressure curve in the barrel and the muzzle velocity of the projectile, values were obtained that were compared with the results of numerical simulations made by Ball-Lagr program.

  4. 4.

    For all simulation variants, i.e., \(7.62 \times 39\) mm FMJ with standard, reduced and increased weights of the powder charge, the values of muzzle velocities, maximum pressures and times to reach maximum pressure are within three standard deviations from the expected values of individual quantities in the distribution T-Student.

A possible further development of the program is to implement other models to solve the main problem of internal ballistics, so that the end user (designer) can choose different models and check which model works best in his conditions. In addition, using the field of optimization techniques, an attempt can be made to implement algorithms that would be able to search for example, the smallest powder weight value at specific powder parameters and a given maximum pressure and projectile muzzle velocity. Additionally, as a result of expanding the program, it would be possible to analyze wave phenomena occurring in the space between the bullet and the cartridge chamber [21,22,23,24,25]. Moreover, the obtained results could be applied as a generalization of the study in the case that the wall of the tube of the gun is considered deformable instead of rigid.